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MULTI-DIMENSIONAL  SIMULATION  OF  ETC  GUN  FLOWFIELDS 

1.0  INTRODUCTION 

This  report  summarizes  progress  made  to  date  in  the  development  of 
advanced  computational  methodology  for  simulating  multi-dimensional, 
electrothermal-chemical  (ETC)  gun  flowfields.  A  research-oriented,  first- 
principles  approach  was  taken  which  utilized  sophisticated,  state-of-the-art 
Roe/TVD  upwind /implicit  numerics,  incorporating  the  relevant  physics  and 
thermochemistry  needed  to  simulate  ETC  processes  in  a  systematic  manner. 
Our  starting  point  for  ETC  simulation  had  involved  the  use  of  the  "original" 
CRAFT  3D  Navier-Stokes  code  [1,2]  whose  applications  had  been  limited  to 
steady-state  jet/propulsive  flows.  CRAFT  was  an  outgrowth  of  the  TUFF  3D 
Navier-Stokes  research  code  developed  by  Molvik  and  Merkle  at  Penn  State 
[3]  under  NASA  Ames  support  for  hypersonic  nonequilibrium  flow  simula¬ 
tion. 

The  National  AeroSpace  Plane  (NASP)  program  prompted  the  devel¬ 
opment  of  a  new  computational  methodology  to  accurately  analyze  chemi¬ 
cally-reacting  flows  with  strong  discontinuities.  TUFF  was  one  of  three  new 
codes  developed  and  the  first  that  was  government-owned.  A  similar  path 
was  taken  by  Rockwell  International  Science  Center  in  their  development  of 
the  "proprietary"  USA  code  [4],  and  somewhat  later  by  Virginia  Pol3d:echnic 
Institute  in  the  development  of  the  GASP  code  [5]  under  NASA  Langley  sup¬ 
port.  The  "1990's"  finite-volume  Roe/TVD  upwind /implicit  numerics  that 
these  codes  were  based  upon  combined  classical  "characteristic  or 
Riemann-invariant  concepts"  with  advanced  matrix  algebra  to  permit  the  ac¬ 
curate  "capturing"  of  strong  discontinuities  (shocks,  contact  surfaces,  flame 
fronts)  in  a  non-oscillatory  manner,  without  the  need  for  stabilizing  artifi¬ 
cial-dissipation  terms.  This  "revolutionary"  breakthrough  in  computational 
fluid  dynamics  was  accompanied  by  the  additional  finding  that  this  class  of 
numerics  was  also  optimal  for  the  treatment  of  "acoustic-driven"  problems 
as  exhibited  by  the  research  of  Beddini  and  students  at  the  University  of 
Illinois  in  analyzing  combustion  instabilities  [6]. 

The  ETC  gun  flowfield  problem  involves  varied  disciplines  as  schema¬ 
tized  in  Figure  1.1  which  shows  the  ingredients  that  the  CRAFT  code  re¬ 
quired  for  simulation.  From  a  numerical  viewpoint,  requisite  baseline  in¬ 
gredients  were  already  in  place.  The  Roe/TVD  formulation  that  CRAFT 
already  contained,  provided  the  ability  to  both  capture  strong  discontinuities 
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Figure  1.1.  ETC  gua  discipHimes  embodied  isa  CRAFT  Navier-Stokes  code, 
in  a  non-oscillatoiy  manner,  and,  to  produce  solutions  of  near-acoustic  accu¬ 
racy  with  negligible  numerical  dissipation.  Earlier  numerical  schemes  (FCT, 
MacCormack)  did  not  provide  this  requisite  behavior.  A  second  principal  in¬ 
gredient  is  the  "fully"  implicit  numerics  —  convective /diffusive/ source  terms 
and  boundary  conditions  are  all  treated  implicitly  which  provides  full-cou¬ 
pling  of  equations  throughout  the  flow  domain.  A  question  raised  in  the 
course  of  this  work  was  the  requirement  to  use  a  fully-implicit  procedure  for 
the  analysis  of  time-accurate  flows  where  large  Courant  numbers  could  de¬ 
grade  the  solution.  The  answer  resided  in  the  need  to  adequately  resolve  all 
flow  regions  and  thus  use  a  non-uniform  grid  where  cell  volumes  could  vary 
by  orders  of  magnitude.  An  implicit  code  can  operate  at  an  "average" 
Courant  number  of  unity,  minimizing  wave  dispersion  errors.  An  explicit 
code  is  restricted  to  operate  at  a  Courant  number  of  unity  for  the  smallest 
cells.  All  other  cells  will  operate  well  below  unity  and  hence,  the  averaged 
Courant  number  can  be  significantly  less  than  one  which  degrades  the 
solution. 
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The  inclusion  of  the  requisite  physical,  thermod^Tiamic,  multiphase 
and  combustion  modeling  parameters  into  CPiAFT  for  ETC  simulation  have 
involved  very  significant  extensions,  as  have  the  specialized  numerical 
upgrades  required  to  simulate  the  d5niamic  interior  ballistic  environment.  It 
has  been  opportune  to  have  had  synergistic  programs  with  several  other 
government  agencies.  These  have  provided  a  framework  for  "piggy-backing" 
the  varied  developments  achieved  in  applying  CRAFT  to  several  different 
problem  areas.  In  addition  to  CRAFT  developmental  support  from  the  U.S. 
Army  Research  Laboratory  for  ETC  gun  simulation  (and  most  recently,  for 
LPG  investigations  -  Ref.  7) ,  additional  support  was  provided  by: 

•  U.S.  Army  Missile  Command  (MICOM),  Redstone  Arsenal,  AL, 
for  tactical  missile  aero /propulsive  simulation  [8-20]. 

•  Naval  Surface  Warfare  Center  (NSWC),  Dahlgren,  VA,  for  mis¬ 
sile  vertical  launcher  system  simulation  [21-24]. 

•  Air  Force  Office  of  Scientific  Research  (AFOSR)/Phillips  Lab- 
Armament  Directorate,  Eglin  AFB,  FL,  for  ram  accelerator 
simulation  [25-28]. 

•  NASA  Langley  Research  Center,  Hampton,  VA,  for  jet/acous¬ 
tics  research  [29-36]. 

•  Air  Force  Wright  Laboratoiy,  Wright- Patters  on  AFB,  OH,  for 
simulation  of  aircraft  plume /wake  interactions  for  signatures 
and  countermeasures  [15,  37-38]. 

The  specific  extensions  to  CFIAFT  in  the  areas  of  numerics,  thermochem¬ 
istry,  multiphase  flow,  and  turbulence  over  the  past  5  years  is  summarized 
below  in  Tables  I  to  IV,  respectively,  including  the  principal  funding  agency 
for  each  extension.  The  CRAFT  EH'C  developmental  work  has  entailed 
problem-specific  extensions  to  CRAFT,  such  as  direct-coupling  with  a 
plasma  capillary  model,  and  varied  extensions  involving  shared  technology 
such  as  large  eddy  simulation  (LES)  turbulence  modeling  for  liquid 
propellants  (shared  with  jet/ acoustic  simulation  activities)  and  Lagrangian 
particulate  modeling  for  solid  propellants  (shared  with  solid  rocket  motor 
simulation  activities). 

In  the  course  of  developing  varied  CRAFT  extensions,  ETC  concepts 
were  changing  rapidly  —  often,  faster  than  developments  could  keep  pace 
with.  Our  starting  point  in  ETC  simulation  involved  analyzing  a  liquid 
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Table  I.  CRAFT  Extensions  /  Numerics 


Table  IV.  CRAFT  Extensions  /  Turbulence 


propellant  working  fluid  with  end  injection  of  the  plasma.  This  changed  to 
analysis  of  central  ullage  tube  concepts  with  "burst-start,"  and  piccolo  tube 
variants.  Liquid  propellant  extensions  focussed  on  the  inclusion  of 
gas /liquid  interactions  and  combustion  modeling  capabilities  into  CRAFT, 
and,  on  simulating  the  large-scale  turbulent  structure  for  which  LES 
methodology  was  found  to  be  extremely  promising. 

As  our  work  with  liquid  propellants  started  to  show  great  promise, 
propellant  concepts  changed  and  we  had  to  shift  gears  quite  rapidly  to 
extend  CRAFT  to  deal  with  solid  propellant  ETC  guns.  Early  work  entailed  a 
simplified  fixed-bed  representation  of  the  solid  propellant.  This  was 
extended  to  a  sophisticated  fluidized  bed  representation  using  a  Lagrangian 
formulation  which  treated  the  interactions  (fluid/thermal)  of  each  deterred 
propellant  ball  discretely.  The  new  formulation  in  CRAFT  was  shown  to 
accurately  simulate  solid  propellant  interior  ballistic  flows  via  detailed  unit 
problem  comparative  studies  with  XKTC  described  in  Ref.  39. 

Our  work  to  date  in  simulating  ETC  gun  flowfields  with  CRAFT  has 
been  described  in  a  number  of  JANNAF  Combustion  Subcommittee  Meeting 
papers  and  related  workshop  presentations  [40-47].  This  report  will  pro¬ 
vide  a  unified  description  of  this  work,  presenting  more  global  details  of  the 
numerical  methodology  not  heretofore  documented.  An  overview  of  CRAFT 
as  a  general  purpose  interior  ballistics  flowfield  solver  with  applications  to 
ETC,  LPG  and  RAMAC  problems;  to  venting/muzzle  blast;  and,  to  projectile 
flyout  and  thermal  heating,  is  described  in  Ref.  48. 
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2,0  CRAFT  CODE  OVERVIEW  AND  SELECTED  STUDIES 

2,1  OVERVIEW  OF  METHODOLOGY  AND  FEATURES 

Interior  ballistic  flowfields  share  common  problems  of  transient 
combustion  with  complex  wave  processes.  The  CRAFT  code  has  been 
applied  to  simulate  interior  ballistic  flowfields  encompassing  electrothermal 
(ETC)  guns,  liquid  propellant  (LP)  guns,  and  ram  accelerators  (as  well  as  to 
simulate  other  propulsive-oriented  problems  as  discussed  earlier).  The 
present  overall  features  of  the  CRAFT  code  are  listed  below  in  Table  V. 


Table  V,  CRAFT  Code  Features 

NUMERICS 

«  1D/2D/AXI/3D  Finite-Volume  Discretization 

®  Implicit,  Higher-Order  Upwind  (Roe/TVD)  Formulation 

«  Fully  Implicit  Source  Terms/Boundary  Conditions 

GRID  FEATURES 

®  Grid  Dynamics  to  Account  for  Moving  Boundaries 

«  Grid  Patching/Blanking  for  Complex  Geometries 

®  Solution-Adaptive  Gridding 

«  Hybrid  Structured/Unstructured  Formulation  for  Multi-Body  Problems 

<»  Real  Gas  Mixtures  (Caloricaily  and  Thermally  imperfect/JANNAF  Thermo 
Tables/Virial  EOS) 

®  Finite-Rate  Chemistry/Arbitrary  Number  of  Species  and  Reactions 

<»  Fully  Implicit  Source  Term  Linearization 

«  Nonequilibrium  Particle/Droplet  Solvers  (Eulerian  and  Lagrangian  Formulations) 

*  Gas/Liquid  Equilibrium  Formulation 

®  Grain/Abtative  Coupling  Including  Surface  Recession  and  Interior  Burning 

o  Coupling  with  3D  Transient  Heat  Conduction  Solution 

®  Generalized  Mass  Transfer  Boundary  Conditions  and  Phase-Change 

®  Coupling  with  3D  Structural  Solver  Cm  Progress) 

TURBULENCE 

«  k-e  Formulation  with  Compressibility/Vortical  Upgrades  and  Several  Low  Re  Near- 
Wall  Formulations 

«  LES  Subgrid  Scale  Models  of  Menon  and  Madabhushi 

•  Particle  Dispersion  Formulations 

APPLICATIONS 

i  «  Electrothermal  Chemical  (ETC)  Gun  -  Uquid  and  Solid 

®  Regenerative  Uquid  Propellant  Gun  (RLPG) 

«  Ram  Accelerator 

<*  Solid/Liquid  Propellant  Rocket  Motors/Exhausts 

®  Ducted  Rocket 

«  Vertical  Launcher  Interactions 

<»  Turbulence/Multi-Phase  Jet  Research 

CFIAFT  is  structured  in  a  finite  volume  framework  and  utilizes  im¬ 
plicit/upwind  numerics  which  entails  sophisticated  and  complex  matrix 
multiplications.  While  transient  interior  ballistic  flowfield  problems  have 
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traditionally  used  explicit  numerics,  the  viewpoint  taken  here  is  that 
implicit  numerics  are  requisite  to  analyze  problems  with  highly  non-uniform 
geometric  scales.  Turbulent  scales  in  near-wall  regions  require  grids  whose 
cell  volumes  are  orders  of  magnitude  smaller  than  average  cell  volumes. 
Explicit  solvers  are  restricted  by  stability  limitations  to  operate  at  time  steps 
associated  with  the  smallest  cell  volume.  Hence,  wave  processes  will  suffer 
from  substantive  dispersive  errors  via  operating  at  very  small  Courant 
numbers.  Implicit  numerics  avoid  such  limitations  and  provide  additional 
stability  for  complex  combusting  flows.  The  Riemann  based  upwind 
numerical  flux  computation  procedure  is  chosen  since  it  permits  the 
accurate  representation  of  strong  gradients  and  discontinuities  (e.g.  shocks, 
flame  fronts,  propellant  grain-ullage  boundaries)  by  aligning  the  numerical 
stencil  with  wave /convective  directions  via  matrix  manipulations. 
Consequently  these  discontinuities  are  captured  sharply  without  the 
characteristic  oscillations  associated  with  earlier  schemes  which  typically 
require  the  use  of  artificial  dissipation  terms  to  suppress  these  spurious 
oscillations. 

The  thermochemistry  in  the  CRAFT  code  has  been  enhanced  from  the 
original  gas-phase  capability  to  include  gas-liquid  mixtures  with  a  discrete 
phase  of  either  droplets  or  propellant  grains.  The  modelling  of  complex 
thermochemistry  associated  with  mixtures  of  gas,  bulk  liquid  and  liquid/ 
solid  particulates  within  this  upwind  numerical  framework  has  required 
significant  developmental  work  to  the  code  and  will  be  described  in  the 
sections  that  follow.  A  number  of  enhanced  grid  capabilities  such  as  grid 
dynamics,  grid  blanking,  and  grid  embedding  have  been  incorporated  to 
better  handle  complex  geometries  whose  volumes  change  and  deforms  in 
time.  The  treatment  of  turbulence  in  CRAFT  is  a  major  philosophical 
departure  from  previous  interior  ballistic  codes.  For  the  highly  transient 
flowfields  of  interest  here,  we  simulate  turbulence  using  a  large-eddy  scale 
(LES)  approach  where  large  vortical  structures  are  captured  as  part  of  the 
unsteady  flow  computation  and  only  the  subscale  turbulence  effects  are 
modeled.  We  note  that  the  development  of  the  code  has  proceeded  in  a 
systematic  fashion  wherein  each  additional  upgrade  has  been  incorporated 
as  necessitated,  guided  by  experimental  data  and  emphasizing  an  improved 
understanding  of  the  physical  phenomena.  In  the  sections  to  follow,  we  will 
describe  the  above  features  and  methodology  in  detail.  In  the  next  section. 
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we  will  describe  the  basic  numerical  framework  and  upwind  flux  algorithm 
as  designed  originally  for  ideal  gas  mixtures.  Thermochemical  upgrades 
required  to  allow  for  treatment  of  non-ideal  gases  and  gas-liquid  mixtures 
within  the  original  upwind  framework  are  described  in  the  subsequent 
section. 

2,2  FINITE-VOLUME  FRAMEWORK  AND  ORIGINAL  GAS  PHASE 
EQUATIONS  FOR  MULTI-SPECIES,  COMBUSTING  FLOWS 

The  conservation  equations  for  multi-species  gas  flows  may  be  written 

in  an  integral  form  for  an  arbitrary  control-volume  as  follows; 

—  \QdV  +  \h^FdS=\DdV  (2,2,1) 

dt  V  ^  y 

where  V  is  the  volume  of  the  control  volume,  n®dS  is  a  vector  element  of  the 
control  surface  with  outward  normal  n  ,  F  represents  both  the  inviscid  and 
viscous  flux  of  the  conserved  quantities  Q  through  the  eontrol  surfaces,  and 
D  contains  any  source  terms.  In  a  finite  volume  framework,  the  arbitrary 
control-volume  is  replaced  with  a  generalized  six  sided  ceil  as  shown  in 
Figure  2.2.1  for  both  a  Cartesian  and  polar  coordinate  system.  For  this  gen- 


Fig,  2,2,1,  Three-dimensiousil  finite-volume  cell, 
eralized  cell  volume,  the  flux  quantity  in  Eq.  (2.2.1)  is  rewritten  as  an  inte¬ 
gral  over  each  face  thereby  yielding 
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///  II  ( 

'^ll(^)+]/2  “  '^y-1/2)^^^^''’  II (^/t+1/2  “ 

(2.2.2) 

~  ^11  (^/+l/2  “  ^i-m ) ^ II  [^j+l/2  ~  ^j-l/2 ) 

Here,  ^  is  the  generalized  streamwise  coordinate,  r\  is  the  normal  coordi¬ 
nate,  and  C  is  the  meridional  coordinate.  The  indices  i,j,k  represent  the  cell 
location.  The  vector  Q  contains  the  conservation  variables  which  are 
defined  at  the  center  of  the  cell  i,j,k.  The  fluxes  are  defined  at  the  cell 
interfaces  which  are  represented  by  non-whole  indices.  The  vectors  E,F, 
and  G  represent  the  inviscid  flux  through  the  cell  interfaces  in  the  tl,  and 
^  directions,  respectively.  R,  S,  and  T  are  the  respective  viscous 
stress/transport  vectors,  and  D  contains  source  terms  resulting  from 
combustion.  We  note  that  since  the  finite-volume  formulation  works  with 
the  physical  cell  surface  and  volume,  we  are  able  to  solve  for  ID,  2D/AXI,  or 
3D  flows  within  a  single  unified  code  by  supplying  the  appropriate  control 
volumes.  Hence,  even  for  a  2D/AXI  flow,  two  grid  planes  in  the  azimuthal 
direction  are  defined  thereby  providing  a  physical  volume  to  the  grid  cells. 

For  a  gas  mixture  of  multi-component  species,  the  conserved  variable 
Q  and  the  inviscid  fluxes  are  defined  as  follows: 


9 


The  metric  quantities  /  .are  components  of  the  cell  face  normal 

L  pointing  in  the  q  direction  with  magnitude  equal  to  the  cell  face  area. 
Similarly,  the  other  metric  quantities  are  the  components  of  the  vectors  M 
and  N,  which  are  the  normal  vectors  to  the  other  families  of  the  cell  faces. 
The  quantity  U  denotes  the  volume  of  the  flux  through  the  cell  face  in  the  ^ 
direction  and  is  defined  as  U  =  £^u  +  £yV  +  (:.w,  while  V  and  W  are  the  corre¬ 
sponding  volume  fluxes  through  the  r\  and  ^  faces.  The  first  five  equations 
solve  for  the  conservation  of  the  mixture  mass  (p),  momentum,  and  total  en- 
ergy  per  unit  volume  (e).  The  subsequent  equations  solve  for  the  partial 
densities  of  (n-1)  chemical  species  (different  species  that  are  accounted 
for).  The  n^h  species  is  partial  density,  pn,  determined  by  the  following 
mass  balance  for  the  mixture  density: 

(2.2.4) 

j=l 

The  equation  of  state  for  the  gaseous  mixture  follows  Dalton's  law  of 
partial  pressures  and  is  written  as  follows: 

(2.2.5) 
M 


Here,  M  is  the  mixture  molecular  weight  which  is  obtained  from  the  indi¬ 
vidual  species  molecular  weights  as  given  below 


M  = 


n 

I 


v^= 


iM  , 

1  s  J 


(2.2.6) 


Here,  Cs  is  the  mass  fraction  of  the  stl^  species  (=ps/p). 
The  expression  for  the  total  internal  energy  is  given  by 


e^p 


h 


-f-  V  +  w 


-p 


(2.2.7) 


where  the  enthalpy  of  the  mixture,  h,  is  determined  by  summing  the  en¬ 
thalpies  of  the  individual  species  as  follows: 


(2.2.8) 


S-\ 

Note  that  the  individual  species  enthalpy  includes  the  heat  of  formation, 
in  addition  to  the  sensible  enthalpy,  hg.  Consequently,  the  energy 
conservation  equation  does  not  contain  a  separate  heat  release  source  term 
since  this  is  implicitly  defined  through  the  heat  of  formation  for  each 
species.  This  point  will  be  further  clarified  in  the  next  section,  when  we 
discuss  the  thermochemistiy  and  temperature  decode  procedure.  We  note, 
as  Eq.  (2.2.5)  indicates,  that  CRAFT  was  originally  configured  to  model 
thermally-perfect  gas  mixtures  only.  The  extensions  to  include  thermally- 
imperfect  virial  equations  of  state  and  a  gas-liquid  equilibrium  formulation 
will  be  described  in  a  later  section. 

2.3  GAS-PHASE  THERMOCHEMISTRY  AND  TEMPERATURE  DECODING 

The  thermochemistry  within  the  CRAFT  code  allows  for  analyzing  a 
mixture  of  calorically  imperfect  gas  species.  Hence,  the  sensible  enthalpy 
and  specific  heat  of  each  species  are  allowed  to  vary  as  functions  of 
temperature  as  denoted  below 


h  =h[T) 
dh 

Cn  =— ^=Cp  (7) 
P,S  JJ’  P.s 


(2.3.1) 


In  the  code,  the  enthalpy  and  specific  heat  are  specified  in  tabular  form  as 
functions  of  temperature.  These  tabular  values  are  then  fitted  with  a  cubic 
spline  to  obtain  a  generalized  functional  relationship  for  the  thermochem¬ 
istry  over  the  entire  temperature  range. 

To  complete  the  thermochemical  specification,  we  still  have  to  de¬ 
code  the  temperature  of  the  gas  mixture  for  given  values  of  internal  energy 
and  species  concentrations.  Since,  the  enthalpy  of  each  species  is  itself  a 
function  of  temperature,  the  internal  energy  becomes  an  implicit  function  of 
temperature.  Therefore,  we  have  to  resort  to  an  iterative  procedure  to  de¬ 
code  the  temperature  from  the  thermodynamic  internal  energy  (e)  and  the 
mass  fractions  of  the  species.  A  Newton-Raphson  procedure  was  adapted 
which  is  illustrated  below. 
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Here,  the  index  k  corresponds  to  the  Newton-Raphson  iterative  index.  The 
iterations  are  continued  until  convergence  is  achieved  i.e.  the  thermody¬ 
namic  internal  energy  is  matched  for  a  given  temperature  value.  We  note 
that  for  a  combusting  case,  the  differences  in  the  heat  of  formation  between 
the  reactants  and  products  ensures  the  correct  heat  release  and  tempera¬ 
ture  increase,  while  ensuring  that  the  internal  energy  is  conserved  exactly. 


2.4  RELEVANT  MATRIX  ALGEBRA 

As  a  precursor  to  our  discussion  of  the  upwind  flux  methodology  and 
the  implicit  solution  procedure,  we  will  introduce  some  of  the  relevant  ma¬ 
trix  manipulations  and  their  physical  interpretations  in  this  section.  To 
simplify  the  algebra,  we  will  deal  with  the  one-dimensional  inviscid  subset  of 
Eq.  (2.2.2)  given  below  in  its  differential  form. 


dQ  dE 


=  0 


(2.4.1) 


While  Eq.  (2.4.1)  does  not  contain  the  viscous  terms  or  the  chemical  source 
terms,  it  is  still  adequate  for  the  purpose  of  mathematically  describing  the 
wave  characteristics  of  the  system.  We  begin  by  defining  the  Jacobian  of  the 
flux  vector  E  as  follows: 

=  (2.4.2) 


The  Jacobian,  A,  can  also  be  viewed  as  a  transformation  matrix  to  linearize 
the  non-linear  flux  vector,  E,  in  time  as  per  Eq.  (2.4.3)  below 


£n+l  =  + 


dQ 


At 


dQ  dt 

.  ^n+\  _  £n  ^  -  2" 


1  2 


(2.4.3) 


Using  Eq.  (2.4.2),  we  rewTite  Eq.  (2.4.1)  in  its  non-consen^ative  form: 


^  + 
dt 


A 


=  0 


(2.4.4) 


Note  that  in  Eq.  (2.4.4),  the  mass,  momentum  and  energy  equations  remain 
strongly  coupled  to  each  other  since  the  Jacobian  A  is  a  full  matrix. 
Therefore,  to  obtain  the  individual  wave  systems  in  the  full  system  of 
equations  we  have  to  decouple  the  components  in  Eq.  (2.4.4). 

The  system  of  equations  in  Eq.  (2.4.4)  are  decoupled  via  the  right  and 
left  eigenvectors  of  the  Jacobian,  A.  We  define  a  matrix,  L,  which  contains 
the  left  eigenvectors  of  the  Jacobian,  A,  while  its  inverse  contains  the  right 
eigenvectors,  R.  By  definition, 

LAR  =  K  (2.4.5) 


where  the  matrix  A  contains  the  eigenvalues  of  the  Jacobian,  A.  To  decouple 
the  equation  system  in  Eq.  (2.4.4)  we  premultiply  it  by  the  matrix  L 


L—+LARL^  =  0 
dt 


Defining  a  new  vector  dQ  =  LdQ,  Eq.  (2.4.6)  may  be  written  as 


(2.4.6) 


dt 


-I- 


(2.4.7) 


In  Eq.  (2.4.7),  the  individual  equations  of  the  system  are  now  decoupled 
from  each  other  since  A  is  a  diagonal  matrix.  Physically,  each  individual 
equation  is  now  a  mathematical  representation  of  a  propagating  wave  whose 
propagation  speed  is  its  eigenvalue  and  the  quantity  that  is  conserved  by  the 
wave  is  the  Riemann  invariant  Q. 

For  the  full  system  of  equations  in  Eq.  (2.2.3)  the  Jacobian,  A,  for  the 
flux  vector  E  is  given  below 
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Similar  expressions  follow  for  Matrices  B  and  C  where  B  =  dFldQ  and 
C=dG/dQ.  These  matrices  are  written  in  a  generalized  form 
accommodating  the  analysis  of  arbitrary  thermodynamics.  The  algebra 
involved  in  determining  matrix  elements  entails  making  derivatives  of  the 
pressure  with  respect  to  various  conservation  quantities  e.g.,  1^1^.  etc.  This 
is  where  the  mixture  thermodynamics  comes  in.  Since  the  matrices  are 
cast  in  generalized  form,  the  inclusion  of  complex  thermodynamics  is 
straightforward  and  is  performed  "off-line."  The  details  of  this  derivation  for 
a  mixture  of  calorically  perfect  gases  are  given  in  the  paper  by  Molvik  and 
Merkle  [3].  The  corresponding  derivation  for  a  gas-liquid  mixture  will  be 
described  in  later  subsection  of  this  section. 

The  eigenvalues  of  the  Jacobian  defined  in  Eq.  (2.4.8)  are  represented 
as 

A  =  DIAG.  {U  +  C,U  -C,U,...,U) 

U  =  iu  +  £v^lw  (2.4.9) 

X  y  z 

c=c^fijnj+jj 

V  X  y  ^ 


where  c  is  the  sound  speed  of  the  gaseous  mixture  and  is  defined  as 


The  eigenvalues  in  Eq.  (2.4.9)  are  associated  with  acoustic  waves  traveling  in 
opposite  directions  propagating  pressure  information  at  the  speed  of  sound, 
and  with  the  waves  convecting  at  the  speed  of  the  fluid  transporting  scalar 
quantities  such  as  entropy  and  chemical  species.  The  elements  of  the  right 
and  left  eigenvectors  for  these  eigenvalues  are  described  in  Ref.  3. 

2.5  UPWIND  FLUX  COMPUTATION 

The  inviscid  fluxes  the  cell  interfaces  are  calculated  via  an  approxi¬ 
mate  Riemann  procedure  described  by  Roe  [49].  In  this  approach,  the  flux 
is  computed  as  the  sum  of  the  contributions  of  individual  waves  crossing  the 
boundary.  The  contribution  of  each  individual  wave  takes  its  propagation 
direction  into  account  thus  providing  the  upwind  characteristics  to  the 
scheme.  The  flux  at  a  cell  interface,  i+1/2,  in  the  ^  direction  (see  Figure 
2.2.1)  is  given  below 

^i+vi  =  +  ^i+i  +  ^+\n.  ~  ^t+\n\ 

Here,  AE'*’  and  AE’  represent  the  flux  changes  associated  with  waves  travel¬ 
ing  in  the  positive  and  negative  directions  respectively.  The  direction  of 
travel  for  each  wave  is  obtained  from  the  eigenvalues  as  discussed  earlier. 
The  flux  differences  AE"^  and  AE‘  are  defined  as 

where  the  matrix  A±|A|  serves  as  an  effective  filter  to  distinguish  between 
positive  and  negative  traveling  waves. 

A  higher  order  inviscid  flux  is  obtained  by  adding  corrective  terms  to 
the  first  order  flux  in  Eq.  (2.5.1).  Since  higher  order  corrections  can  induce 
spurious  oscillations,  these  terms  are  checked  to  ensure  that  they  satisfy  the 
Total  Variation  Diminishing  (TVD)  criteria  [50].  The  higher  order  flux  is 
written  as 
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(2,5.3) 


where  (t)=l/3  for  a  third  order  scheme.  The  characteristic  variable  differ¬ 
ence,  (Aa),  is  the  quantity  which  is  checked  for  TVD  criteria.  The  charac¬ 
teristic  variable  is  defined  as 


'^“(•+1/2  "  h+yi  ( 2/+1  ) 


(2.5.4) 


The  flux  limiters  are  defined  using  minmod  operators  as  follows: 

(o=)  =minmod[(®).^j^,P(o).^3^] 


(2.5.5) 


where  the  minmod  operator  is  defined  as 

minmod[j:,  v]  =  sign(.r)*max[0,min{|4>'*sign(j:)}]  (2,5.6) 

and  P  is  the  compression  parameter  that  is  restricted  to  lie  in  the  range 

l<P<2zi  (2.5.7) 

l-(t) 


For  additional  details  on  TVD  schemes  the  reader  is  referred  to  Ref.  50. 


2.6  IMPLICIT  SOLUTION  PROCEDURE 

To  formulate  the  implicit  solution  procedure,  we  begin  by  linearizing 
the  inviscid  flux  terms,  the  viscous  terms,  and  the  source  terms  with  re¬ 
spect  to  the  conserved  variables  Q.  The  linearized  terms,  when  combined, 
yield  a  sparse  block-pentadiagonal  matrix  for  2D  cases  and  a  block-septadi- 
agonal  matrix  for  3D  problems.  These  sparse  matrices  are  not  inverted  ex- 
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actly  since  it  is  computationally  expensive,  but  instead  an  approximate 
Alternate-Direction-Implicit  (ADI)  procedure  is  employed.  The  linearization 
procedure  for  the  inviscid  flux  terms  and  the  resulting  ADI  operator  are  de¬ 
scribed  below. 

The  first  order  linear  flux  in  Eq.  (2.5.1)  is  rewritten  as 


/+1/2  9  L  ^+1  i 


where 
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(2.6.2) 


Linearizing  Eq.  (2.6.1)  by  making  use  of  Eq.  (2.4.3)  we  obtain 
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where  A  0  is  the  change  in  the  conserved  variable  and  is  defined  as 

A!2,=(e"^^-!2")  (2.6.4) 


The  flux  terms  in  the  q  and  ^  directions  can  similarly  be  linearized  to  yield 
expressions  corresponding  to  Eq.  (2.6.1).  The  viscous  terms  and  the 
chemical  source  terms  can  be  linearized  in  a  similar  fashion  as  well.  For  the 
sake  of  notational  simplicity,  we  will  denote  the  linearized  Jacobians  for 
these  terms  as  W  without  going  into  the  details  which  are  provided  in  Ref.  3. 

Combining  the  linearized  terms  we  obtain 
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Ae^j=-A(i 


dQ  dE  dF  dG  ,, 

- j - j_ - j - Y  —  £) 

dt  aq  ac  tTi.c 


(2.6.5) 


Here  the  operator  5^Ai  is  defined  as 

yBiji,  =  +  yx  -  <2.6.6) 


while  the  corresponding  operators  5r|Bj  and  5^Ck  are  defined  as 


(2.6.7) 


^  {^M/2  ~  "  ^k-\/2^^i.J.k-i 


As  Eq.  (2.6.5)  indicates,  the  block  matrix  on  the  left  hand  side  is  a  septadi- 
agonal  (since  we  are  formulating  for  the  3D  equation  system).  To  invert  this 
large  matrix  in  a  cost  effective  manner,  we  employ  the  ADI  procedure. 

The  ADI  operator  inverts  the  large  septadiagonal  matrix  (or 
pentadiagonal  matrix  in  a  2D  case)  by  splitting  the  operator  as  a  product  of 
three  independent  sweeps  in  the  p,  and  ^  directions.  The  ADI  approxi¬ 
mation  for  the  implicit  matrix  is  given  as 
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We  note  that  the  operator  in  each  sweep  is  now  only  a  tridiagonal  matrix 
which  can  be  inverted  efficiently  using  the  Thomas  algorithm.  Furthermore, 
while  the  ADI  operator  introduces  an  error  into  the  solution  from  the  ap¬ 
proximations  it  entails,  this  error  can  be  reduced  by  performing  iterations  at 
each  time  step  and  converging  to  a  time-accurate  solution  [51]. 


2.7  TURBULEMCE  MODELING:  kc  AND  LES 

Turbulent  flows  are  characterized  by  the  existence  of  spatially  and 
temporally  fluctuating  fields  with  a  wide  range  of  length  and  time  scales. 
Depending  on  how  the  various  length  and  time  scales  are  resolved,  three 
approaches  are  currently  used  in  the  numerical  simulation  of  turbulent 
flows:  the  time-averaged  or  Re3molds-Averaged  (RANS)  approach,  the  Direct 
Numerical  Simulation  (DNS)  approach,  and  the  Large  Eddy  Simulation  (LES) 
approach.  In  the  RANS  approach,  the  ensemble-averaged  Navier-Stokes 
equations  are  solved  along  with  turbulence  closure  models  for  the  Reynolds- 
stresses.  The  RANS  turbulence  models  are  calibrated  using  experimental 
observations  and  often  incorporate  heuristic  corrections.  RANS  is  the 
approach  utilized  for  simulating  most  steady  or  slowly  varying  flows.  The 
large  scale  features  of  flowfields  differ  considerably  from  one  flow  to 
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another,  e.g.,  ducted  flows,  aerodynamic  flows,  jet/wake  flows,  stratified 
flows,  etc.,  all  have  significantly  different  flow  structures.  RANS  models  use 
a  single  scale  equation, to  represent  a  broad  range  of  turbulent  length  scales. 
The  invariance  of  modeling  coefficients  cannot  be  made  very  general  for 
differing  flows  with  any  reliability.  For  flows  with  rapid  temporal  variations 
such  at  ETC,  use  of  RANS  time-averaged  assumpbons  are  very  questionable. 

In  DNS.  the  non-averaged /time-dependent  Navier-Stokes  equations 
are  solved  with  resolution  of  all  temporal  and  spatial  scales  of  the  flow 
required  eliminating  the  need  for  turbulence  modeling.  The  results  from 
the  DNS  approach  are  very  accurate  and  reveal  valuable  information  on  the 
turbulence  structures.  Since  the  direct  numerical  simulations  resolve  all 
length  scales,  the  computational  demands  are  very  large  (requiring 
hundreds  of  hours  of  CPU  time  on  super  computers).  Also,  because  of  the 
scaling  of  CPU  time  and  memory  with  powers  of  Reynolds  number,  direction 
numerical  simulations  are  limited  to  analyzing  low  Reynolds  number  flows. 

A  compromise  to  the  DNS  approach  is  the  technique  of  Large  Eddy 
Simulations  (LES)  in  which  the  large  scales  of  turbulence  are  directly 
simulated  while  the  small  (inner)  scales(subgrid  scales  or  SGS)  are  modeled. 
The  large  scales  represent  the  anisotropic  part  of  the  energy  spectrum  and 
contain  most  of  the  energy.  The  small  scales  are  more  isotropic  in  nature 
and  relatively  independent  from  the  resolved  part  of  the  spectrum  (grid 
scales).  The  results  from  the  LES  can  be  quite  accurate  because  the  only 
approximation  involved  is  the  modeling  of  the  SGS  terms  which  do  not 
contain  much  energy.  LES  has  the  advantage  of  reduced  empiricism  over 
the  RANS  approach  and  is  applicable  to  flows  with  rapid  temporal  variations. 
It  does  not  require  the  fine  grids  entailed  in  DNS  simulation  and  is  thus 
applicable  to  higher  Reynolds  number  flows. 

2.7.1  ke  Model 

A  popular  turbulence  model  in  the  RANS  framework  is  the  two- 
equation  k£  model  in  which  partial  differential  equations  are  solved  for  the 
transport  of  turbulent  kinetic  energy  (k)  and  turbulent  dissipation  (e),  and 
the  eddy  viscosity  is  related  to  these  two  quantities.  The  basic  high 
Reynolds  number  form  of  the  ke  equations  [Launder  et  al.,  1972]  (Ref.  52)  is 

given  by: 
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where  the  turbulent  produetion,  P,  is  given  by: 
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(2.7. 1.2) 


and  the  turbulent  stress,  pujUj,  is  evaluated  using  the  eddy  viscosity  assump¬ 
tion: 
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The  turbulent  viscosity,  pt.  is  defined  by: 


(2.7. 1.4) 


and  the  "standard"  coefficients  are  as  follows: 

C^  =  09,  C.  =  L45,  C2=t9,  a^=tO.  a^  =  t3 

The  presence  of  solid  boundaries  affects  turbulence  characteristics  in  the 
near  vicinity  of  the  wall.  The  basic  high  Reynolds  form  of  the  kc  equation 
has  to  be  modified  in  order  to  extend  its  applicability  to  the  wall.  This 
involves  the  addition  of  "low  Reynolds  number"  terms  to  the  equations 
whose  influence  is  strongest  close  to  the  wall  and  diminishes  away  from  the 
wall.  In  the  CRAFT  eode,  the  low  Reynolds  number  variant  of  Chien  [53]  has 
been  incorporated.  In  this  model,  the  source  terms  are  modified  as: 

k  eqn:  P-pe — (2,7. 1.5) 
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where 
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and  the  coefficients  are  now  given  by: 


(2.7. 1.7) 

(2.7. 1.8) 


q=1.35,  C2=1.8,  C3=.0115.  ^=0.5,  C^=0.09 

Here  L  is  the  distance  from  the  closest  wall  and  |i  is  the  laminar  viscosity. 
The  ke  equations  are  solved  in  a  coupled  manner  with  the  rest  of  the 
equations  (k  and  ke  are  added  to  the  Q  array  and  matrix  sizes  are  expanded). 
The  ke  model  in  CRAFT  has  been  used  primarily  for  steady  or  slowly 
changing  flow  calculations. 


2.7.2  Large  Eddy  Simulation  (LES) 

For  high  Re3molds  number  flows  of  practical  interest,  direct  numerical 
simulation  resolving  all  spatial  and  temporal  scales  is  not  feasible  and  the 
turbulence  must  be  modeled.  Conventional  time-averaged  Reynolds-stress 
(RS)  modeling  for  unsteady  flows  requires  that  turbulent  time-scales  are 
small  in  comparison  to  convective  time  scales.  Such  averaging  may  be  ade¬ 
quate  for  slowly  accelerating  flows  or  flows  with  low  frequency  periodic  be¬ 
havior.  It  is  not  generally  adequate  for  short-duration  transient  flows  and 
flows  with  higher  frequency  periodic  behavior.  For  such  flows,  the  time 
scales  of  the  large  scale  turbulent  structure  and  the  convective  processes 
are  comparable.  A  large-eddy  simulation  (LES)  procedure  provides  a 
workable  method  for  transient  flows  whereby  the  large  scale  turbulent 
structure  is  directly  simulated  by  the  flow  solver  and  the  finer  scale, 
dissipative  structure  (which  cannot  be  resolved  numerically  with  practical 
grids)  is  modeled.  The  modeling  involves  a  spatial  filtering  process  which 
involves  subgrid  stress  (SGS)  terms  as  will  be  discussed  below.  Table  VI 
summarizes  the  basic  distinctions  between  RANS  and  LES  modeling  ap¬ 
proaches. 
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LES  vs  RANS  OVERVIEW  /  INCOMPRESSIBLE 


1.  DECOMPOSITION 
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3.  TERMS  FROM  AVERAGING 


o 

LES 

um".  -f- 

UM  .  -  UM  .  H-  U.U".  ~\-UM  . 

INVARIANT? 

^  J 

[  i  J  i  J  J  [  i  J  1  J  J 

o  RANS  u'.u'.  ^  NOT  INVARIANT,  COEFFICIENTS  ARE 

'  J 

PROBLEM-DEPENDENT 


Reliable  SGS  models  for  non-homogeneous  turbulent  compressible  flows 
are  still  in  the  developmental  stage  [54,55].  Issues  such  as  the  backscatter 
of  turbulent  kinetic  energy  from  small  scales  to  the  large  scales,  etc.,  need 
further  investigation  especially  when  high  Re3molds  number  flows  are 
simulated  with  moderate  resolution.  In  studies  described  in  this  report,  a 
simple  subgrid  model  which  is  a  compressible  extension  of  the  Smagorinsky 
model  [56]  is  used  where  the  effect  of  subgrid  scales  is  assumed  to  be  only 
dissipative. 

The  simulation  of  turbulence  using  LES  methodology  begins  with  the  fil¬ 
tering  of  the  small  scale  effects  from  the  large  scale  motion  in  the  full 
Navier-Stokes  equations.  The  filtering  for  Favre- averaged  variables  is  given 
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as 


Here,  G.ix.-X.  .A.)  is  the  filter  function  and  is  linked  to  the  filter  width,  Ai. 
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The  time  varying  variable  d  is  now  expressed  as 


<|)  =  0  +  (j)' 


(2.7.2.2) 


where  ^  is  the  Favre -averaged  large  scale  component  of  (p,  and  0',  is  the 
small  scale  component. 

Appl5dng  the  filtering  described  above,  the  Navier-Stokes  equations  are 
expressed  as 
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The  terms  superscripted  SGS  denote  the  small  scale  terms  which  have  to  be 
modelled. 

The  SGS  terms  in  the  momentum  equation,  when  expanded,  can  be 
viewed  as  a  sum  of  three  terms:  Leonard  Stress,  Cross  Stress,  and  Reynolds 
Subgrid  Stress  as  given  below 
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Here,  L,  C,  and  R,  are  the  Leonard,  Cross  (forward  scattering)  and  Reynolds 
subgrid  (backward  scattering)  terms,  respectively.  These  three  terms  are 
modeled  using  a  simplified  compressible  generalization  of  the  Smagorinsky 
model  [56]  which  assumes  that  the  subgrid  scale  turbulence  kinetic  energy 
is  smaller  than  the  thermod3mamic  pressure.  The  simplification  yields 
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where  Sy  is  the  strain  tensor  of  the  Favre -averaged  large  scale  component, 
given  by 
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and  the  subgrid  eddy  viscosity  is  given  by 
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Here,  A  is  the  effective  filter  width,  given  by 
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(2.7.2.81 


A^.Ay.A^  are  the  filter  widths  in  the  x,  y,  and  z  directions. 

In  our  use  of  CRAFT,  the  filter  widths  are  approximated  by  the 
respective  grid  spacings.  A  value  of  0.02  is  used  for  Cr.  The  subgrid  terms 
in  the  energy  equation  are  given  by 
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The  SGS  terms  in  the  energy  equation  are  modeled  using  scale  similarity 
yielding 
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where  the  subgrid  kinetic  energy  is  neglected  when  accounting  for  the 
fluctuations  of  the  total  energy. 


By  definition,  filtering  is  implicit  in  any  numerical  simulation  using  a 
finite  number  of  grid  points  because  length  scales  smaller  than  the  grid  size 
cannot  be  resolved.  This  implicit  filtering  is  called  a  top-hat  or  box  filter  in 
finite-difference  representation  and  a  sharp  cut-off  filter  in  a  spectral 
(Fourier)  representation.  In  addition  to  the  implicit  filtering,  an  explicit 
filter  can  also  be  applied,  such  as  the  popular  Gaussian  filter,  with  filter 
widths  different  from  the  local  grid  spacing.  This  explicit  filtering  is  usually 
referred  to  as  "prefiltering"  and  the  argument  in  its  favor  is  that  since  the 
length  scale  used  in  the  SGS  model  is  decided  by  the  filter  width,  it  would 
be  prudent  to  make  the  SGS  model  independent  of  the  numerical  resolution 
by  applying  a  prefilter. 

2.8  SPECIALIZED  GRID  FEATURES 

In  this  section  we  highlight  upgrades  to  the  grid  methodology  in 
CRAFT  which  enhances  its  capabilities  for  geometrically  complex,  dynamic 
volume,  interior  ballistic  flowfields.  The  three  primary  upgrades  are  as 
follows: 

a)  Dynamic  Grids 

b)  Grid  Blanking 

c)  Grid  Embedding 

In  the  following  paragraphs,  we  will  briefly  describe  the  salient  features  of 
each  of  these  three  upgrades. 

2.8.1  D3mamic  Grids 

Djoiamic  grid  techniques  allow  for  the  motion  and  expansion  of  the 
cell  volume  as  the  flow  domain  changes  in  time.  The  grid  movement  is 
coupled  to  the  flow  solution  and  the  additional  fluxes  as  well  as  the 
volumetric  effects  generated  from  the  movement  of  the  grid  are 
incorporated  into  the  numerical  procedure.  Our  methodology  in  CRAFT 
closely  follows  the  formulation  of  Vinokur  (Ref.  57)  and  we  refer  the  reader 
to  his  report  for  details. 

The  modified  form  of  the  integral  conservation  equations  for  a  time- 
varying  grid  are  given  as 
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(2.8.1. 1) 


We  note  that  in  comparison  with  Eq.  (2.2.1),  Equation  (2. 8. 1.1)  reveals  an. 
additional  integral  in  time  since  the  cell  volumes  and  areas  are  changing. 
The  discretized  form  of  Eq.  (2.8. 1.1)  after  taking  the  additional  time  integral 
into  account  is  written  as  follows  (Ref.  57): 
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The  modifications  arising  from  the  dynamic  grid  movement  are  classi¬ 
fied  into  three  sets  of  terms.  Term  I  generates  a  volumetric  source  term  for 
the  conserved  variables.  This  term  ensures  that  provided  there  is  no  net 
flux  or  combustion,  the  change  in  volume  proportionally  alters  the  mass  and 
energy  per  unit  volume.  Term  II  contains  the  changes  to  the  flux  vectors  (E, 
F,  and  G)  which  now  have  to  account  for  additional  flux  being  generated  by 
the  movement  of  each  face  of  the  generalized  finite-volume  cell.  For  in¬ 
stance,  the  flux  vectors  E  and  F  are  written  as 
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Here  Ur  is  the  additional  volumetric  flux  generated  by  the  movement  of  the 
^  face,  while  Vj-  is  the  corresponding  movement  of  the  rj  face.  These  addi¬ 
tional  fluxes  due  to  face  movement  are  computed  using  the  following  rela¬ 
tionship: 
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Here  SV^  is  the  volume  swept  by  the  c  face  over  the  period  of  the  time  step. 
Similar  relationships  are  also  defined  for  the  r\  and  ^  faces  as  well.  The 
metrics  that  are  used  to  evaluate  the  flux  vectors  etc.)  are  the  aver¬ 

age  surface  area  of  the  cell  over  the  time  step.  Finally,  Term  111  contains  the 
volumetric  changes  to  the  source  terms. 

To  evaluate  the  upwind  flux  terms  (Ei+1/2)  at  the  cell  interfaces,  we 
recompute  the  eigenvalues  of  the  modified  flux  vectors.  The  new  eigenval¬ 
ues  are  derived  to  be 


A  =  D1AG.{U'  +  C,U'-  C, U\... U')  (2.8. 1.5) 

where 

U'  =  U-Uj.  (2.8. 1.6) 

Hence  by  modifying  the  definition  of  the  volume  flux  through  each  face,  the 
upwind  flux  relationships  derived  in  Section  2.5  may  be  used  unchanged. 
The  d5niamic  grid  formulation  described  above  has  been  tested  extensively 
to  ensure  that  the  code's  conservation  capabilities  do  not  deteriorate  when 
strong  grid  movement  and  distortion  are  present. 

2.8.2  Grid  Blanking 

Grid  blanking  increases  the  versatility  of  the  CRAFT  code  for  problems 
involving  complex  geometries.  The  grid  blanking  feature  in  CRAFT  has  been 
adapted  from  the  PARC  code  [58]  and  works  in  conjunction  with  the  ADI 
procedure  to  invert  the  matrix  arrays  in  the  implicit  sweep.  Grid  blanking 
is  achieved  by  using  a  patching  algorithm  which  logiCcdly  decomposes  a  grid 
containing  internal  boundaries  into  a  family  of  patches.  In  addition  to  sim¬ 
plifying  gridding  issues,  grid  blanking  also  permits  a  very  generalized 
boundary  specification  procedure  since  each  boundary  plane  can  now  be 
broken  up  into  segments  having  different  boundary  conditions. 

The  grid  blanking  feature  in  CRAFT  is  illustrated  by  a  computation  of 
the  venting/ muzzle  blast  from  a  gun  barrel  (see  Figure  2.8.2. 1).  One  end  of 
the  gun  barrel  is  shut  while  the  other  end  is  opened  to  the  atmosphere  at 
t=0.  The  transient  flow  inside  the  barrel,  and  the  external  plume  are 
computed  concurrently  using  grid  blanking.  The  single  rectangular  grid 
implemented  for  the  problem  is  shown  in  Figure  2. 8. 2. 2.  It  is  clear  that 
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since  the  gun  barrel  lies  in  the  interior  of  the  grid,  analyzing  this  problem 
requires  the  use  of  grid  patching  since  two  sets  of  boundary  conditions  are 
required  -  one  for  the  interior  surface  of  the  barrel  wall,  and  one  for  the 
exterior  surface  of  the  barrel  wall  which  is  exposed  to  the  ambience.  The 
grid  patches  used  for  this  problem  are  shown  in  Figure  2. 8. 2. 3. 
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Fig.  2. 8.2.3.  Grid  patches  generated  for  flux  and  boundary  condition 
calculations. 

In  Figure  2. 8. 2. 4,  we  track  the  exhaust  from  the  gun  barrel  by  plotting 
the  stagnation  enthalpy  contours  at  1,  2,  and  3  ms  respectively.  The  initial 
condition  of  the  gas  in  the  barrel  has  a  pressure  and  temperature  ratio  10 
times  that  of  the  respective  ambient  conditions  outside.  At  1  ms,  the  high 
pressure  gas  in  the  barrel  sends  an  underexpanded  plume  out  creating  a 
blast  wave.  By  2  ms,  the  plume  expands  further  out  thereby  reducing  the 
pressure  within  the  barrel  since  there  is  only  a  finite  amount  of  gas  within 
the  barrel.  By  3  ms,  the  gas  within  the  barrel  overexpands  below  the  pres¬ 
sure  of  the  plume  outside.  The  higher  back  pressure  sends  a  compression 
back  into  the  barrel  and  flow  reversal  occurs  which  is  clearly  seen  in  the 
stagnation  enthalpy  contours.  In  problems  with  finite -thickness  solid 
volumes,  grid  points  within  the  volume  are  blanked  out.  In  our  recent  work 
on  the  ram  accelerators  [27,28]  a  coupled  3D  transient  conduction  capability 
has  been  added  so  that  the  transient  wall  temperature  can  be  calculated  as 
part  of  the  basic  computational  procedure.  The  solid  grid  points  blanked 
out  now  serve  to  solve  the  conduction  equation. 

2.8.3  Grid  Embedding 

The  grid  embedding  feature  is  necessary  for  unsteady  flows  where  the 
flow  domain  is  changing  substantially  in  volume  as  in  interior  ballistic 
problems.  This  feature  is  used  in  conjunction  with  the  dynamic  movement 
of  the  grid.  If  a  grid  cell  expands  and  becomes  larger  than  a  user-specified 
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Fig.  2. 8. 2. 4.  Stagnation  enthalpy  contours  plotted  at  three  time  levels  during  venting 
of  the  olume. 


limit,  it  is  subdivided  into  smaller  cells.  This  ensures  that  as  the  domain  of 
interest  increases,  the  grid  continues  to  resolve  spatial  gradients  accurately. 
The  flow  properties  in  the  divided  cell  are  obtained  by  interpolating  from 
the  original  coarser  grid  each  time  a  division  occurs.  To  illustrate  this 
feature,  we  plot  the  initial  grid  in  a  ballistic  chamber  in  Figure  2.8.3.1a.  In 
Figure  2.8.3.1b,  we  plot  the  corresponding  grid  in  the  chamber  later  in  the 
ballistic  cycle  at  which  point  the  domain  has  expanded  considerably.  As  is 
apparent  even  though  the  spatial  domain  shows  a  three  fold  increase,  the 
ability  to  capture  discontinuities  or  to  resolve  wave  motion  is  unimpaired 
because  the  grid  embedding  procedure  has  added  an  adequate  number  of 
grid  points  to  maintain  the  original  spatial  accuracy. 


Fig.  2.8.3. 1.  Grid  inside  gim  barrel  at  two  time  levels  exhibiting 
grid  embedding  procedure. 

2.9  BOUNDARY  CONDITIONS 

In  this  section,  we  discuss  the  boundary  condition  procedures  in 
CRAFT  that  are  applied  in  a  manner  consistent  with  the  finite  volume  formu¬ 
lation.  We  focus  attention  here  on  inviscid  flow  conditions  since  viscous 
boundary  layers  were  generally  not  resolved  in  our  ETC  gun  simulations. 
However,  generalized  viscous  wall  conditions  are  operational  in  CRAFT  and 
we  refer  the  reader  to  our  ram  accelerator  publications  [28],  where  we 
describe  the  analysis  of  wall  viscous  effects  for  transient  flows. 


The  inviscid  boundar>'  condition  methodolog}'  is  based  on  the  method 
of  characteristics.  In  this  procedure,  we  determine  the  number  of  waves 
entering  the  flow  domain  from  the  outside  at  each  boundary,  and  apply 
boundary  conditions  to  account  for  the  information  they  convect  in.  For  the 
remainder  of  the  waves  which  reach  the  boundary  from  within  the  flowfield, 
the  information  is  computed  from  the  equations  of  motion.  The  direction  of 
propagation  for  each  wave  is  deteirnined  from  its  eigenvalue  i.e.  a  positive 
eigenvalue  denotes  a  wave  traveling  in  the  positive  direction  and  vice  versa. 

The  classes  of  inviscid  boundary  conditions  that  are  of  interest  in  CTC 
interior  ballistic  flow  problems  include  subsonic/supersonic  inflows  and  wall 
reflection /centerline  conditions.  The  boundary  condition  for  supersonic  in¬ 
flows  is  straightforward  since  all  the  eigenvalues  are  positive  and  hence  all 
the  characteristics  (or  waves)  are  entering  the  boundary  from  outside  the 
flow  domain.  Therefore,  all  flow  variables  are  specified  as  boundary 
conditions  and  the  flux  at  the  boundary  face  is  completely  defined.  Subsonic 
inflow  boundary  conditions  are  more  involved,  since  one  eigenvalue  (U  -  C) 
is  negative  while  the  remaining  four  (U+C,  U,  U,  U)  are  positive.  Hence  four 
boundary  conditions  are  specified.  Typically  the  following  four  conditions 
are  applied  for  steady  inflows: 

( 1 )  Total  Enthalpy  (Ht) 

(2)  Total  Pressure  (Pt) 

(3,4)  Cross-flow  angularity  (V/U  and  W/U)  which  is  not  typically 
obtainable  and  assumed  to  be  zero  in  the  plane  parallel  to  the  in¬ 
flow  boundary  condition  (V,  W)  which  ensures  that  the  inflow  is 
perpendicular  to  the  boundary  plane. 


The  fifth  condition  which  corresponds  to  the  negative  characteristic  is 
obtained  by  extrapolating  the  characteristic  variable  corresponding  to  this 
wave.  This  is  accomplished  through  the  solution  of  the  folloAving  equation 
between  the  boundary  and  the  interior  point: 


dQ  dE 
dt  ^ 


=  0 


(2.9.1) 


Here  the  row  of  the  left  eigenvector  matrix  L  which  corresponds  to  the 

(U  -  C)  eigenvalue.  The  boundary  conditions  result  in  five  nonlinear 
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equations  which  are  solved  in  a  coupled  matrix  form  via  a  Newton-Raphson 
iterative  procedure  to  obtain  the  flow  quantities. 

For  unsteady  applications,  the  two  conditions  on  Pt  and  Ht  are  not 
utilized  since  these  are  appropriate  only  at  steady  state.  These  two 
conditions  are  replaced  with 

T^  =  0,  L  —  =  0  (2.9.2) 

^  ^  dt 


Equation(2.9.2)  allows  for  waves  to  pass  through  the  inflow  and  provides  a 
non-reflective  boundary  condition. 

Inviscid  wall  boundary  as  well  as  centerline  symmetry  conditions  are 
handled  by  generating  ficbtious  cells  outside  of  the  computational  domain  so 
the  cell  interface  flux  that  coincides  with  the  wall  or  the  s)rmmetry  plane 
can  be  handled  as  if  it  were  an  interior  plane.  The  flow  variables  of  the  cells 
inside  the  domain  are  reflected  about  the  symmetry  plane  to  produce  a  mir¬ 
ror  cell  on  the  other  side.  The  scalar  flow  variables  in  the  reflected  cell  (e.g. 
density  and  energy)  are  identical  to  the  interior  cell  values  and  are  given  by 


p''  =  p 

e’^  =  e 


(2.9.3) 


On  the  other  hand,  the  gas  velocities  in  the  two  cells  are  mirror  images  to 
ensure  that  the  net  normal  mass  flux  through  the  boundary  is  zero  although 
there  still  is  a  pressure  flux.  The  velocities  in  the  reflected  cell  are  given  by 

u’^  =  u-2i  ii  u  +  i  v  +  i  w) 
x\  X  y  z  } 

v^  =  v-2t  ii  u  +  i  v  +  i  w]  (2.9.4) 

y\  X  y  z  ) 

w’^  =  w  —  2i  (i  u  +  i  v  +  i  w] 

X  >>  z  J 


The  unit  metrics  in  the  above  expressions  Eire  those  of  the  boundary  cell 
under  consideration.  This  transformation  preserves  the  magnitude  of  the 
velocity  while  reorienting  its  direction. 


2.10  SELECTED  FUNDAMENTAL  STUDIES 

In  this  section,  we  describe  several  fundamental  studies  which 
demonstrate  the  ability  of  the  CRAFT  code  to  accurately  analyze  problems 
involving  finite-rate  combustion,  turbulence  with  large  scale  vortical 
structures,  and,  transient  wave  processes  in  interior  ballisbc  chambers.  We 
begin  by  describing  some  very  basic  studies  for  shock  tubes  with  high 
pressure  ratios  to  emphasize  the  ability  of  CRAFT  to  capture  strong  shocks 
without  spurious  oscillations.  In  the  next  subsection,  we  describe 
fundamental  combustion  studies  in  premixed  hydrogen-air  mixtures  where 
we  focus  on  the  need  to  accurately  resolve  the  spatial  structure.  This  is 
followed  by  simulations  of  classical  shear  layer  experiments  with  large 
density  gradients  where  large  scale  turbulent  vortical  structures  were  first 
identified  experimentally.  Finally,  we  discuss  some  preliminary  simulations 
of  ETC  like  configurations  where  high  pressure  gas  is  injected  into  a  gun 
chamber  and  the  resulting  wave  proeesses  are  tracked. 

2.10.1  Shock  Tube  Studies 

The  ability  of  CRAFT  to  simulate  unsteady  flowfield  phenomena  is 
assessed  by  computing  one-dimensional  solutions  for  a  shock  tube  sealed  at 
both  ends.  A  shock  tube  10  units  in  length  was  selected  with  a  diaphragm 
located  midway,  separating  the  high-pressure  chamber  on  the  left  from  the 
low-pressure  chamber  on  the  right.  Initial  studies  were  conducted  for  a 
pressure  ratio  of  10  across  the  diaphragm.  Figure  2.10.1.1  shows  a  time 
history  of  the  pressure  distribution  along  the  shock  tube.  Figure  2.10.1.2 
shows  a  corresponding  plot  for  a  subsequent  study  conducted  with  an  initial 
pressure  ratio  of  1000.  In  both  cases,  the  code  was  found  to  operate  very 
robustly  with  minimal  numerical  diffusion,  and  a  shock  pattern  which  did 
not  attenuate  with  time.  Further,  it  was  found  that  the  code  operated  very 
satisfactorily  at  large  CFL  numbers  -  a  feature  which  is  directly  attributable 
to  the  implicit,  strongly-coupled,  higher-order,  upAvind  numerical  formula¬ 
tion  of  CRAFT. 


2.10.2  Gas-Phase  CombnistioEi  Studies 

For  applications  (such  as  the  Ram  Accelerator),  where  the  principal 
phenomena  is  premixed  gas-phase  combustion,  the  accuracy  of  the  numer¬ 
ical  simulation  hinges  on  the  code's  ability  to  capture  flame  fronts  as  well  as 
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shock  waves.  We  have  performed  some  fundamental  combustion  studies  to 
evaluate  the  ability  of  the  CRAFT  code  to  capture  flame  fronts  in  premixed 
combusting  flows  using  as  test  cases,  the  shock  induced  combustion 
experiments  of  Lehr  [59].  Two  sets  of  calculations  have  been  performed. 
The  first  for  a  steady  combustion  front  and  the  second  for  an  oscillating 
flame  front.  Recently,  these  cases  have  been  numerically  simulated  by  other 
authors  [60,61]  using  different  algorithms  as  well  as  different  procedures  for 
evaluating  the  chemistry.  CRAFT  utilizes  strongly-coupled  methodology  with 
a  "full"  linearization  of  the  chemical  producbon  term,  a.,  with  respect  to  the 
Q  array.  Computing  these  test  cases  provides  a  rigorous  evaluation  of  the 
ability  of  CFiAFT  to  compute  such  flows. 

The  first  test  case  simulated  is  for  a  steady  shock-induced  combustion 
front.  The  geometry  is  a  15mm  diameter  sphere-cylinder,  moving  at  Mach 
6.46  into  a  stoichiometric  mixture  of  hydrogen  and  air.  The  premixed  gas  is 
at  a  pressure  of  0.42  atm  and  a  temperature  of  292  K.  The  simulation  was 
first  performed  on  a  relatively  coarse,  uniform  grid  of  71  x  71  as  shown  in 
Figure  2.10.2.1a.  The  reaction  systems  utilized  are  summarized  in  Table  Vll 
[62,63].  The  computed  density  contours  for  the  uniform  grid  are  plotted 


a.  Uniform  Grid  h.  Density  Contours 


Fig.  2.10.2.1.  IJTffilfoTrm  nMmerical  grid  and  computed  density  contours 
shown  for  steady  shock-induced  combustion. 


plotted  in  Figure  2.10.2.1b  while  the  experimental  values  for  the  shock  and 
flame  front  are  plotted  as  circles.  As  is  e\adent  from  Figure  2.10.2.1b,  with  a 
uniform  grid  the  numerical  calculation  underpredicts  the  induction  zone 
significantly,  away  from  the  centerline.  A  similar  observation  has  also  been 
made  by  Sussman  [61].  Both  sets  of  reactions  were  found  to  produce 
eomparable  results. 

_ Table  Vn.  HYDROGEN/AIR  RATES _ 
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To  improve  upon  the  computed  predictions  without  increasing  the 
grid  size,  the  uniform  71  x  71  grid  was  adapted  using  SAGE  [64]  emplo3dng 
density  as  the  controlling  variable.  The  final  adapted  grid  is  shown  in  Figure 
2.10.2.2a  and  the  computed  results  using  this  grid  are  shown  in  Figure 
2.10.2.2b.  The  results  for  the  adapted  grid  compare  well  with  the  experi¬ 
mental  values  -  the  induction  zone  spacing  is  captured  accurately  as  it 


a.  Adapted  Grid 


b.  Density  Contours 


Fig.  2.10,2.2.  Adapted  grid  and  computed  density  contours  shown  for 
steady  shock-induced  combustion. 

broadens  in  the  streamwise  direction.  This  improvement  is  attributed  solely 
to  the  additional  grid  points  in  the  induction  zone  resulting  from  adaption 
and,  again,  the  results  with  both  reaction  sets  were  comparable.  Sussman 
[61]  states  that  in  addition  to  grid  adaptation,  he  had  to  modify  the  treat¬ 
ment  of  the  chemical  source  term  to  accurately  compute  the  induction  zone. 
Our  experience  with  CRAFT,  indicates  that  when  an  appropriately  adapted 
grid  is  used,  the  induction  zone  can  be  captured  without  any  modifications. 
Note  the  very  sharp  capturing  of  the  flame  front  by  CRAFT. 

The  second  test  case  simulates  an  unsteady,  shock-induced  combus¬ 
tion  flowfield.  The  geometry  and  mixture  are  the  same  as  in  the  previously 
described  steady  case,  but  the  freestream  Mach  Number  now  has  a  lower 
value  of  4.79.  In  the  experiment,  high  frequency  oscillations  are  observed  in 
the  shadowgraphs  (Figure  2.10.2.3).  This  oscillatory  behavior  is  thought  to 
be  the  result  of  combustion  instabilities  caused  by  the  interaction  between 
the  bow  shock  and  the  compression  waves  generated  by  the  energy  release. 
A  detailed  discussion  of  the  proposed  instability  mechanisms  can  be  found  in 
Ref.  60. 


Fig.  2.10.2.3.  Shadowgraph  of  unsteady  M=4.79  case  from  Lehr. 

The  simulation  was  performed  on  an  evenly  spaced  166  x  121  grid 
utilizing  the  first  reaction  system[62]  employed  in  the  steady  case.  Second- 
order  time-accuracy  was  employed  for  this  unsteady  calculation.  Figure 
2.10.2.4a  shows  density  contours  at  one  point  in  time  illustrating  the 
"pulsating"  behavior  of  the  energy  release  front.  Figure  2.10.2.4b  presents 
an  X  vs  t  diagram  of  density  contours  plotted  along  the  centerline  between 
the  bow  shock  and  the  projectile  nose.  This  figure  illustrates  the  periodic 
movement  of  the  energy  release  zone.  The  estimated  frequency  of  this 
oscillation  is  approximately  450-500  KHz  which  is  lower  than  the 
experimental  value  of  720.  These  results  are  similar  to  those  obtained  by 
Wilson  &  Sussman  [60].  They  demonstrated  the  sensitivity  of  the  oscillation 
frequency  to  the  H-O2  chain  branching  reaction  and  were  able  to  improve 
upon  the  predicted  frequency  by  utilizing  the  rate  expressions  of 
Jachimowski  [63].  We  plan  to  repeat  this  unsteady  calculation  with  these 
rates  to  assess  this  sensitivity. 


Brown  &  Roshko[65]  were  amongst  the  earliest  researchers  to 
identify  large  scale  coherent  structures  in  shear  layers  with  density 
gradients  and  their  data  has  formed  the  basis  for  many  unsteady  simulations 
(see,  e.g.,  the  Euler  simulation  of  Chien  [66]  with  excitation  forced  at  the 
dominant  instability  modes).  The  case  of  interest  is  described  in  Figure 
2.10.3.1  (for  helium  and  nitrogen  streams)  which  shows  the  time-averaged 


Fig.  2.10.3.1.  Unsteady  shear  layer  simulation  with  large  density 

variations:  a)  schematic;  and  b)  time-averaged  solution. 

structure  of  the  shear  layer  in  a  reduced  transverse  domain  (the 
computational  domain  extends  well  beyond  what  is  displayed).  The 
calculation  is  continuously  processed  to  ascertain  when  time-invariance  is 
achieved  (via  averaging  over  -1,000  steps)  and  at  what  axial  position 
asymptotic  similarity  is  achieved.  Conventional  Favre-averaging  is  imple¬ 
mented  to  obtain  mean  and  rms  turbulent  variables /stresses  as  summarized 
below  in  Table  VIII  which  omits  the  subscale  stress  terms. 


Table  Vin.  Density-Weighted  (Favre)  Averaging  of  Unsteady  Flow  Solutions 


SOLUTION  =>  Q  (Xj,  t)  =  [p,  pu,  pv,  pw,  pe|.l 

PRIMITIVE  VARIABLES 
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Figure  2.10.3.2  shows  instantaneous  contours  of  nitrogen  mass 
fraction  and  clearly  indicates  the  presence  of  unsteady  large  scale  turbulent 
structures.  Referring  to  Figure  2.10.3.1,  a  stable  asymptotic  state  is 
achieved  at  the  vertical  position  indicated.  Figure  2.10.3.3  exhibits  the  time 
variation  of  axial  velocity  at  4  points  in  this  shear  layer  (see  Figure  2.10.3.1 
for  the  positions).  Figure  2.10.3.4  exhibits  the  Favre  and  Reynolds  time- 
averaged  mean  axial  velocity  profile  at  this  axial  position  and  at  several 
positions  downstream  (the  downstream  solutions  still  exhibit  some 
unsteadiness  indicating  that  the  calculation  is  not  yet  fully  time-invariant). 
Figure  2.10.3.5  exhibits  axial  rms  velocity  fluctuations  profiles  (ir). 
Analogous  results  for  instantaneous  density  fluctuations  at  the  same  4  points, 
and,  for  the  rms  density  {'p')  fluctuation  profiles  at  the  same  4  axial  stations 
are  shown  in  Figures  2.10.3.6  and  2.10.3.7,  respectively.  The  qualitative 
characteristics  of  the  results  obtained  are  quite  satisfying  and  were  obtained 
with  low  frequency/low  amplitude  forcing  with  and  without  the  utilization  of 
a  subscale  stress  model  (differences  were  found  to  be  negligible). 
Quantitative  comparisons  with  time-averaged  density  fluctuation  data  are 
quite  good  as  shown  in  Figure  2.10.3.8. 
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2.10.4  Gas/Gas  Gun  Chamber  Studies 

Simulations  of  high-pressure,  plasma  injection  into  a  gun  chamber 
and  subsequent  projectile  acceleration  were  performed  at  early  stages  of  our 
ETC  work  as  an  initial  assessment  of  numerics.  For  the  case  described,  both 
the  plasma  and  propellant  were  represented  by  air  and  a  non-combusting 
firing  was  simulated.  The  plasma  was  injected  into  the  gun  mixing  chamber 
and  a  total  duration  of  ~2ms  was  simulated.  During  the  plasma  injection 
process,  a  complex  flow  structure  was  observed  in  the  chamber  with  strong, 
embedded  discontinuities.  Multiple  shock  reflections  occurred  and  a  Mach 
disc  formed  downstream  of  the  plasma  injector.  An  unsteady  shear  layer  is 
found  to  emanate  from  the  Mach  disc  triple  point  and  separates  the  hot  core 
from  the  tube  walls.  The  plasma  shock  wave  travels  to  the  base  of  the  pro¬ 
jectile  protrusion  and  minor  reflections  are  observed.  However,  most  of  the 
blast  wave  travels  along  the  projectile  protrusion  unbl  it  reflects  from  the  far 
end.  This  initiates  projectile  motion  while  the  reflected  shock  rapidly  trav¬ 
els  back  to  interact  with  the  Mach  disc  structure.  At  this  point  (the  end  of 
simulation),  the  injector  is  close  to  "unstarting."  Contours  of  pressure  and 
temperature  are  shown  in  Figures  2.10.4.1  and  2.10.4.2  at  various  time  in¬ 
stances  to  provide  snapshots  of  the  flowfield  while  Figure  2.10.4.3  shows  the 
variation  of  pressure  along  the  centerline  of  the  mixing  chamber  at  different 
times. 
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Newer  gun  design  concepts,  such  as  the  Electrothermal-Chemical 
(ETC)  gun  and  the  Regenerative  Liquid  Propellant  Gun  (RLPG),  have  utilized 
liquid  propellants,  either  wholly  or  in  conjunction  with  solid  propellants,  to 
provide  enhanced  performance.  The  attraction  of  liquid  propellants  lies  in 
their  high  chemieal  energy  content.  Furthermore,  they  offer  a  potential  for 
better  control  and  tailoring  of  the  pressure-time  curve  since  the  interaction 
between  the  gas  and  liquid  phases  plays  a  dominant  role  in  the  combustion 
process.  For  example,  one  of  the  proposed  designs  for  the  electrothermal- 
chemical  gun  involves  using  bulk  liquid  propellant  whose  combustion  is  con¬ 
trolled  by  injecting  high  pressure  plasma  which  ignites  and  energizes  the 
liquid  propellant.  The  transient  mixing  processes  at  the  interface  of  the 
liquid  and  gaseous  plasma  are  very  complex  and  involve  the  generation  of 
droplets  from  the  bulk  liquid  phase.  The  simulation  of  such  complex  multi¬ 
phase  processes  requires  an  integrated  code  which  permits  a  gas-phase,  a 
bulk-liquid  phase,  and  a  discrete  droplet  phase  to  coexist,  with  each  phase 
potentially  being  in  non-equilibrium  with  the  other  phases.  A  first-princi¬ 
ples  simulation  of  droplet  formation  in  this  environment  is  possible  at  a  fun¬ 
damental  level,  but  would  entail  substantive  research  and  very  specialized 
solution  adaptive  gridding  to  make  it  into  a  practicable  component  of  an 
engineering-oriented  code.  Empirical  models  of  droplet  formation  [67]  are 
not  catered  to  the  ETC  environment  and  no  relevant  modeling  work  in  this 
area  has  been  performed  to  date.  In  view  of  the  state-of-the-art  having  such 
inadequacies,  an  approximate  gas-liquid  equilibrium  approach  was  formu¬ 
lated  as  a  first  step  towards  the  subsequent  development  of  a  more  general¬ 
ized  non-equilibrium  approach  (to  be  described  in  Section  6).  The  gas-liq¬ 
uid  equilibrium  formulation  is  described  in  what  follows. 


3,1  GAS-LIgUID  EgUBLIBRIUM  FORMULATION 

The  basic  premise  of  the  gas-liquid  equilibrium  formulation  is  the  as¬ 
sumption  of  dynamic  and  thermal  equilibrium  between  the  two  distinct 
phases  i.e.  the  gas  and  liquid  phases  have  identical  velocity  and  temperature 
at  any  spatial  point  in  which  they  coexist.  The  equilibrium  assumption  al¬ 
lows  us  to  combine  the  conservation  equations  for  each  individual  phase  to 


obtain  a  unified  set  of  equations  for  the  mLxture.  The  equation  system  is 
completed  by  enforcing  Amagat's  law  between  the  phases  -  the  gas  and  liq¬ 
uid  phase  pressures  are  identical  while  the  partial  volumes  .sum  up  to  the 
cell  volume.  The  conservation  equations  are  supplemented  by  the  equations 
of  state  for  the  gas  and  liquid  phases.  The  gas-phase  has  a  generalized  virial 
equation  of  state  to  account  for  compressibility  effects  at  high  pressures,  as 
will  be  discussed  in  detail  later.  The  liquid  equation  of  state  can  be  speci¬ 
fied  either  in  a  tabular  form,  or  as  an  analytical  function.  The  basis  of  the 
gas-liquid  equilibrium  formulation  is  summarized  in  Table  IX.  The  method¬ 
ology  described  has  been  adapted  from  procedures  employed  by  Coffee  [68] 
for  RLPG  simulation. 

Table  IX.  Basis  of  Gas/Liquid  Equilibrium  Formulation 


Pl-Ul.Tl,Vl 
LIQUID _ 

Pg  ’  ^G’  ^G’  ^G 

GAS 

•  TWO-PHASE  MIXTURE  IS  GOVERNED  BY  AMAGATS  LAW 

TOTAL  VOLUME  =  SUM  OF  PARTIAL  VOLUMES 
't'G  +  <1>L  =  1 

PRESSURE  IS  IDENTICAL  IN  BOTH  PHASES 
Pg=Pl 

•  GAS  AND  LIQUID  PHASES  TAKEN  TO  BE  IN  EQUILIBRIUM 

0g=0l 

Tg=Tl 

THIS  ALLOWS  US  TO  FORMULATE  GLOBAL  CONSERVATION  EQUATIONS  FOR  THE  MIXTURE 
RATHER  THAN  SEPARATE  EQUATIONS  FOR  EACH  PHASE 


The  conservation  equations  for  the  gas-liquid  mixture  are  written  in 
vector  form  as 


j^gll(‘^y-i-l/2  ‘^7-1/2)'^^^^ 


(3.1.1a) 


Here,  the  nomenclature  used  is  the  same  as  that  of  the  previous  section  with 
Q  being  the  vector  of  dependent  variables,  while  E,  F,  and  G  are  the  flux 
vectors  in  the  and  r\  directions,  respectively.  These  vectors  are  given  as 
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(3. 1.1b) 


The  metrics  the  volumetric  fluxes  U,  V,  and  W,  and  the  mix¬ 

ture  total  internal  energy,  e^,  have  the  same  definitions  as  in  Section  2.  As 
before,  we  formulate  the  equations  for  n  general  species.  However,  unlike 
the  pure  gas-phase  formulation,  we  now  have  two  phases  present  and  we 
designate  the  nth  species  to  be  the  bulk  liquid  phase.  The  remainder  of  the 
(n-1)  species  are  the  various  gaseous  chemical  species  which  are  individually 
modelled.  The  first  five  equations  are  the  conservation  equations  for  the 
gas /liquid  mixture  where  pm  is  the  mixture  density  and  is  defined  as  follows: 


(3.1.2) 


In  Eq.  (3.1.1b),  (j)g  is  the  volumetric  fracdon  of  the  gas,  while  pg  and  p^  are 
the  gas  and  liquid  densities.  Since  the  gas  phase  is  composed  of  (n-1) 
different  species,  we  define  pg  as 

n-1 

P  =  XP-  (So1o3) 

i-\ 

We  define  a  relative  gas  mass  fraction,  yi,  as  follows: 

.v,=—  = 
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(3.1.4) 


Note  that  the  relative  gas  mass  fractions,  y^,  are  distinct  from  the  species 


mass  fractions 


C. 


The  mixture  enthalpy  is  given  by 


p  h  =p(t)/z+p,  l-d  /j 


n-I 

I-'', 

/  =  1 


hf  +  h.(T) 


(3.1.5) 


The  conservation  equations  in  Eq.  (3.1.2)  are  supplemented  by  the 
equations  of  state  for  the  gas  and  liquid  to  complete  the  system.  The  pres¬ 
sure  of  the  gas  is  the  sum  of  the  individual  gas  species  and  is  given  as 

=  (3.1-6) 

(  =  1 


The  pressure  of  each  individual  gas  species  is  specified  using  the  non-ideal 
(thermally  imperfect)  virial  equation  of  state  as  follows: 


p.RT 

P  .=^ - 

M. 


1  +  50  -h  C  p 


(3.1.7) 


The  coefficients  By  and  Cy  are  the  second  and  third  virial  coefficients  for  the 
gaseous  mixture.  Typically,  the  second  and  third  virial  coefficients  are  avail¬ 
able  for  individual  species.  However,  the  methodology  to  deduce  the  virial 
coefficients  for  a  gaseous  mixture  from  these  individual  values  is  not  well 
defined.  Analytical  expressions  have  been  derived  for  binary  mixtures  of  two 
species:  however,  such  expressions  involve  defining  additional  cross-coeffi¬ 
cients  which  are  difficult  to  obtain.  In  view  of  these  uncertainties,  the  virial 
coefficients  for  a  generalized  mixture  of  n  species  have  been  approximated 
by  a  simplified  mass  weighted  averaging  of  individual  species  coefficients 
(69)  as  given  below 
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(3.1.8) 


The  liquid  equation  of  state  is  specified  either  in  tabular  form  or  as  an 
analytical  function,  if  available.  In  the  studies  to  be  described,  we  have  used 
the  following  analytical  expression  employed  by  Coffee  [68]: 


where  is  the  liquid  density  at  a  reference  state,  k^,  is  the  bulk  modulus  at 
zero  pressure,  and  the  index  k2  is  a  measure  of  the  liquid's  compressibility. 

To  complete  the  thermochemical  specification,  we  have  to  obtain  the 
gas  volumetric  fraction,  ^g.  and  the  temperature  of  the  mixture  from  the  con¬ 
served  variables,  q.  As  in  the  pure  gas  phase  case,  the  decoding  procedure 
is  an  iterative  process.  However,  the  gas /liquid  problem  is  more  compli¬ 
cated  since  we  have  to  simultaneously  iterate  for  two  variables  (T  and  (j)g) 
rather  than  just  the  temperature.  The  two  conditions  iterated  on  are  as 
follows: 

m  p  «  p  2^ 

(3.1.10) 


Po 


(3.1.9) 


The  Newton -Raph son  iterative  procedure  for  Eq.  (3.1.10)  is  written  symbol¬ 
ically  as 
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(3.1.11) 


y-r+i  =  tP  +  {^dTf 

♦/*'=♦/+(<;♦  J' 

Eq.  (3.1.11)  is  iterated  on  until  the  error  on  the  right  hand  side  is  reduced 
to  machine  accuracy,  thereby  5rielding  the  mixture  temperature  and  gas  vol¬ 
umetric  fraction. 
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3.2  MATRIX  MODIFICATIONS 

The  inclusion  of  the  gas-liquid  thermochemistry  described  above  into 
the  upwind /implicit  methodology  is  complex  and  requires  the  re-evaluation 
of  the  flux  Jacobians,  as  well  as  the  eigenvalues  and  eigenvectors.  As  we  had 
previously  discussed  in  Section  2.4,  the  Jacobian  of  the  flux  vectors  requires 
the  evaluation  of  the  derivatives  such  as  ,  etc.  To  evaluate  the  deriva- 

m 

tives  of  pressure  we  begin  by  looking  at  the  functional  dependence  of 
pressure 


Inspection  of  Eq.  (3.2.1)  reveals  that  while  the  pressure  differential  depends 
on  (d(t)g),  (t)g  itself  is  not  a  variable  in  the  vector  Q.  Therefore  the  primary  ob¬ 
stacle  in  determining  the  pressure  derivative  is  in  eliminating  the  depen¬ 
dence  of  the  pressure  differential  on  (d(j)g).  This  is  done  by  enforcing  the 
differential  form  of  Amagat's  law, 


dP  =dP, 
S  L 


(3.2.2a) 


Eq.  (3.2.2),  is  substituted  back  into  Eq.  (3.2.1)  to  obtain  the  pressure  differ¬ 
ential  in  terms  of  the  variables  solved  in  the  vector  Q  as  follows; 


dP  =  dP  =dP, 
s 


HUp^,dp  ,dT 


(3.2.2b) 


The  details  of  the  Jacobian  A  are  given  in  Appendix  A. 

The  acoustic  speed  of  the  two-phase  system  is  determined  from  the 
eigenvalues  of  the  Jacobian  (see  Appendix  B  for  details)  and  is  given  by 
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(3.2.3) 


In  the  above  expression 
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(3.2.5) 


As  Eq.  (3.2.3)  indicates,  the  acoustic  speed  in  a  two-phase  gas/liquid  mix¬ 
ture  behaves  differently  than  in  either  a  pure  gas  or  a  pure  liquid.  For  two- 
phase  mixtures,  the  acoustic  speeds  of  the  individual  phases  do  not  combine 
in  a  linear  fashion,  but  instead  exhibit  a  harmonic  relationship.  The  varia¬ 
tion  of  the  acoustic  speed  with  the  volumetric  fraction  of  the  gas  (for  a  given 
pressure  and  temperature)  has  a  "bath-tub"  shape:  while  the  values  at  the 
two  limits  are  the  respective  single  phase  values,  the  acoustic  speed  drops 
sharply  away  from  these  limits.  The  strong  dependence  of  the  acoustic 
speed  on  the  volumetric  fraction  puts  more  stringent  requirements  on  the 
numerical  algorithm  since  it  is  now  required  to  propagate  waves  with  the 
correct  speed  in  mixtures  where  the  two-phase  composition  is  var3dng 
rapidly  (e.g.  across  the  plasma/liquid  propellant  mixing  layer.  The  upwind 
methodology  in  CRAFT  enables  us  to  compute  these  waves  accurately,  as  we 
shall  demonstrate  for  two-phase  shock  tube  unit  problems  later  in  this  sec¬ 
tion. 

3.3  PLASMA  COUPLIMG  PROCEDURE 

To  simulate  ETC  gun  systems,  we  have  to  model  the  injection  of 
plasma  into  the  chamber.  The  plasma,  which  initiates  and  sustains  the  com¬ 
bustion  of  the  propellant,  is  a  critical  component  of  the  ETTC  gun  system. 
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The  mass  flux  and  properties  of  the  plasma  being  injected  depend  upon  the 
instantaneous  pressure  field  within  the  gun  chamber  since  the  plasma  is 
subsonic  for  a  substantial  part  of  the  injection  cycle.  Hence,  to  compute  ac¬ 
curate  injection  properties  for  the  plasma,  the  plasma  generation  code  must 
be  coupled  to  the  CRAFT  code.  The  plasma  generation  code  used  in  our 
ETC  simulations  has  been  developed  by  Powell  at  ARL  (see  Ref.  70  for 
details).  The  code  for  plasma  generation  is  one-dimensional,  quasi-steady, 
and  isothermal,  and  has  been  found  to  provide  adequate  results  for  ETC 
simulation  for  varied  conditions. 

It  responds  "instantaneously"  to  the  chamber  conditions  at  the  injec¬ 
tor  Interface  and  cannot  deal  with  very  high  pressures  that  would  induce 
"backflow"  (flow  from  chamber  into  plasma  generator).  Under  such  condi¬ 
tions,  injector  "door"  is  shut  temporarily.  The  coupling  procedure  between 
CRAFT  and  Powell's  code  is  illustrated  in  Table  X.  This  coupling  procedure 
has  been  verified  by  comparing  the  Mach  number  at  the  injection  plane  for 
the  coupled  solution  with  that  computed  by  the  plasma  code  given  measured 
back  pressures  from  an  experimental  firing.  Such  a  Mach  number  compari¬ 
son  is  shown  in  Fig.  3.3.1.  These  results  indicate  that  the  coupling  proce¬ 
dure  is  functioning  adequately. 


Table  X.  Plasma  Coupling  Procedure 


VELOCITY,  TEMPERATURE  AND 
PRESSURE  (IF  CHOKED) 
COMPUTED  BY  PLASMA 


I{t) 

CODE/SPECIRED  AS  INJECTION 

PLASMA 

CONDITION  TO  CRAFT 

CRAFT 

CURRENT  SPECIFIED 

AS  INPUT  FROM 

CODE 

U  T  P 

WpLi  'pL»  >^pL 

CODE 

P(t) 


PRESSURE  AT  COUPLING  PLANE 
COMPUTED  BY  CRAFT  SPECIFIED 
AS  BOUNDARY  CONDITION 
TO  PLASMA  CODE 


1.1 


3.4  PROPELLANT  COMBUSTION/DECOMPOSITION  MODEL 

The  combustion/ decomposition  of  bulk  liquid  propellant  in  CRAFT  is 
modelled  via  a  two-step  process.  In  the  first  step,  the  bulk  liquid  is  con¬ 
verted  to  an  intermediate  gaseous  form  which  subsequently  burns  in  the 
gaseous  phase  to  generate  the  final  products.  The  rate  for  both  steps  is 
specified  in  finite-rate  Arrhenius  form.  The  numerical  values  for  the  rate 
coefficients  were  obtained  after  performing  one-dimensional  numerical  ex¬ 
periments  to  determine  an  appropriate  rate  of  pressure  rise.  This  tempera¬ 
ture  dependent  combustion  formulation  has  two  implications  from  a  physical 
viewpoint.  The  first  observation  we  make  is  that  the  mixing  between  the 
hot  plasma  and  the  colder  liquid  controls  the  mixture  temperature  and  con¬ 
sequently  the  overall  rate  for  the  first  step  of  the  decomposition  process 
(e.g.  this  process  is  essentially  diffusion-controlled).  The  second  is  that  the 
temperature  dependent  decomposition  rate  allows  us  to  prevent  the  exis¬ 
tence  of  liquids  at  very  high  temperatures.  The  liquid  is  modeled  to  de¬ 
compose  at  a  finite  rate  when  its  temperature  goes  above  the  boiling  point. 
As  we  discussed  earlier,  the  actual  process  by  which  the  bulk  liquid  com¬ 
busts  is  extremely  complex  and  involves  both  temperature  dependent  de¬ 
composition,  as  well  as  burning  of  droplets  which  are  generated  at  the 
phases  interface.  The  complex  physics  at  this  interface  is  not  well  under¬ 
stood,  and  quantitative  data  to  model  it  from  a  more  fundamental  basis  is  not 
currently  available.  The  approach  we  have  followed  is  a  credible  approxima- 
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tion  that  attempts  to  incorporate  the  dominant  effects  of  mLxing  between 
the  hot  plasma  and  the  bulk  liquid. 

For  the  ETC  simulations  to  be  discussed  later  in  this  section,  a  global 

four-component  (n=4)  formulation  has  been  utilized  comprised  of: 

1 .  Gaseous  Plasma  ( pj ) 

2.  Gaseous  Reactants  (p.,) 

3.  Gaseous  Products  (p^) 

4.  Bulk  Liquid  Propellant  (p^) 

The  plasma  is  inert  and  transfers  energy  to  the  working  fluid  via  convective 
and  conduction  processes  only.  Radiative  transfer  has  been  neglected  (for 
simplicity)  although  its  contributions  are  felt  to  be  significant  by  several  in¬ 
vestigators  (see,  e.g.  the  comments  of  Gillian,  Ref.  71). 

The  finite-rate  description  of  the  decomposition  and  gas-phase  com¬ 
bustion  processes  are  as  follows: 


Decomposition: 

(PROPELLANT) (PROPELLANT)y^p 


K 

V 


=  A  r 

V 


-B  ) 

V 


(3.4.1) 


Combustion: 


30  (PROP) 


10  C6»2  +  14/^2+51 


K  =AT 

c  c 


k-BA 


(3.4.2) 


The  specific  rate  constants  utilized  for  our  ETC  gun  predictions  will  be 
given  in  a  later  subsection. 


3.5  FUNDAMENTAL  STUDIES 

In  this  section,  we  describe  some  fundamental  studies  to  demonstrate 
the  accuracy  and  range  of  application  of  the  gas-liquid  formulation  in  CRAFT. 
We  begin  by  looking  at  shock  propagation  in  high  pressure  (non-ideal)  gases 
by  computing  shock  tubes  where  the  pressure  ratio  is  constant  but  the  mean 
pressure  level  increases.  This  is  followed  by  a  validation  exercise  for  the 
gas-liquid  formulation  in  which  we  compute  a  series  of  two-phase  shock 
tube  cases.  The  results  of  CRAFT  are  compared  with  the  calculations  per¬ 
formed  by  Coffee  [68]  using  a  different  numerical  methodology.  Finally,  to 
highlight  the  versatility  of  the  gas-liquid  formulation,  we  simulate  the  de- 


formation  and  aerodynamic  breakup  of  a  cylindrical  slug  of  bulk  liquid  by 
mixing  with  the  surrounding  air. 


3.5.1  High-Pressure  Gas  Shock  Tube  Study 

To  illustrate  the  effect  of  high  pressure  compressibility  in  non-ideal 
gases  [as  represented  by  the  virial  equation  of  state,  Eq.  (3.1.9)],  numerical 
shock  tube  studies  were  performed  at  varied  pressures.  Figure  3.5. 1.1  ex¬ 
hibits  a  10:1  pressure  shock  tube  solution  for  three  different  pressure  levels 
(10:1  atm,  100:10  atm,  1000:100  atm).  The  shock  tube  was  taken  to  be  of 
unit  length  and  is  plotted  along  the  x-axis  in  Fig.  3. 5. 1.1,  while  the  y-axis  is 
the  time  scale.  The  pressure  plots  in  Fig.  3. 5. 1.1  show  the  temporal  evolu¬ 
tion  of  the  shock  and  expansion  waves  as  they  reflect  off  the  ends  of  the 
shock  tube.  We  note  that  at  the  lowest  pressure  levels  (10:1  atm),  non-ideal 
gas  effects  are  negligible  and  the  wave  propagation  is  similar  to  that  of  an 
ideal  gas.  As  the  mean  pressure  increases,  the  non-ideal  effects  start  be¬ 
coming  important  and  the  acoustic  speeds  as  well  as  shock  propagation 
speeds  go  up.  The  increase  in  the  propagation  speeds  is  evidenced  by  the 
greater  number  of  shock  reflections  for  the  highest  pressure  level  case. 
Without  the  inclusion  of  the  virial  EOS,  all  these  predictions  would  be  the 
same,  resulting  in  very  significant  errors  for  high  pressures. 
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3.5.2  Validation  Studies  for  Two-Phase  Shock  Tubes 


To  validate  the  gas/liquid  equilibrium  formulation  in  CRAFT,  a  two- 
phase  shock  tube  problem  was  chosen  as  a  test  case.  ,  The  solutions 
computed  by  CRAFT  were  compared  with  those  obtained  by  Coffee  using  a 
MacCormack  algorithm  based  code  [68].  Four  cases  were  computed,  all 
having  pressure  ratios  of  2:1,  with  a  high  pressure  of  100  MPa  and  a  low 
pressure  of  50  MPa.  Table  XI  gives  the  details  of  each  case. 

Table  XI.  Gas /Liquid  Test  Cases 


SOLUTIONS  FOR  TWO-PHASE  1-D  SHOCK  TUBES  WERE 
COMPARED  WITH  TERRY  COFFEE'S  COMPUTATIONS 


CASE 

1: 

Liquid:  high  pressure  side 
Gas:  low  pressure  side 

CASE 

II: 

Liquid:  low  pressure  side 
Gas:  high  pressure  side 

CASE 

III: 

Uniform  mixture  of  gas 
and  liquid  on  both  sides 

CASE 

IV: 

Radial  shock  tube  - 
uniform  mixture  with  high 
pressure  at  core 

Figure  3.5.2. 1  compares  instantaneous  snap-shots  of  the  pressure 
profiles  at  0.02  ms  and  0.04  ms  for  Case  1.  In  this  calculation,  the  high 
pressure  side  contains  liquid  and  the  low  pressure  side  contains  gas.  Both 
the  CRAFT  upwind  solution  and  the  MacCormack  central  difference  solution 
(with  no  added  artificial  dissipation)  compute  the  correct  shock  amplitude 
and  propagate  both  the  shock  and  the  rarefaction  at  similar  speeds.  The 
MacCormack  solution,  as  expected,  shows  numerical  oscillations  around  the 
shock  front,  while  the  upwind  formulation,  formulated  to  capture  disconti¬ 
nuities,  provides  a  smooth  solution.  The  corresponding  Case  2  solution  with 
the  gas  on  the  high  pressure  side  and  the  liquid  on  the  low  pressure  side  is 
shown  in  Fig.  3. 5.2.2.  The  pressure  profiles  again  compare  quite  well  and 
provide  confidence  that  the  new  upwind  formulation  is  computing  the  cor¬ 
rect  solution  for  these  two-phase  flowfields. 
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Fig.  3. 5.2.2.  Pressture  profiles  for  a  two-phase  shock  tube:  gas  =  high 
pressrure  side  and  liquid  =  low  pressure  side. 

The  next  series  of  calculations  (Cases  3  and  4)  were  performed  for  a 
uniform  gas /liquid  mixture  where  the  volume  fraction  of  the  gas  is  0.501. 
An  axial  shock  tube  case  is  simulated  in  Case  3,  and  the  corresponding  radial 
case  for  Case  4.  Figure  3. 5. 2. 3  compares  pressure  profiles  at  0.05  and  0.15 
ms  for  the  axial  shock  tube  problem  of  Case  3.  The  shock  amplitudes  agree 
very  well.  The  wave  propagation  speeds  which  determine  the  location  of  the 
shock  show  minor  differences  at  the  larger  time  of  0.15  ms.  This  difference 
was  found  to  be  attributable  to  use  of  a  virial  equation  of  state  (EOS)  in  the 
CRAFT  upwind  formulation,  and  a  Noble-Abel  EOS  with  a  co-volume  term  in 
the  Coffee/MacCormack  procedure.  The  two  EOS's  give  very  similar  results 
at  lower  pressures.  However,  at  larger  pressures,  higher  order  virial  terms 
have  to  be  considered  for  the  virial  EOS  to  match  the  Noble-Abel  EOS.  For 


these  test  cases,  only  the  first  virial  coefficients  were  specified  and  hence 
the  acoustic  wave  speeds  are  somewhat  different  for  the  two  formulations. 
The  pressure  histories  at  the  two  ends  are  shown  in  Fig.  3. 5. 2. 4.  The  ampli¬ 
tude  and  the  frequencies  compare  well  taking  into  account  the  small  differ¬ 
ences  in  the  acoustic  speeds. 


The  pressure  oscillation  comparisons  in  the  Case  4  radial  shock  tube 
are  shown  in  Fig.  3. 5. 2. 5.  As  in  Case  3,  the  fluid  on  both  the  high  and  low 
pressure  sides  is  a  uniform  gas/liquid  mixture.  Unlike  the  axial  case,  the 
pressure  waves  in  the  radial  shock  tube  problem  have  a  strong  amplification 
due  to  axisymmetric  effects  near  the  centerline.  The  pressure  oscillations 
on  both  the  centerline  as  well  as  the  upper  wall  are  exhibited.  The  oscilla¬ 
tions  do  not  attenuate  significantly  in  the  duration  of  the  calcultion.  A  peri¬ 
odic  limit-cycle  with  a  pressure  amplitude  of  30-40  MPa  is  generated  at 
larger  times. 
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Fig.  3.5.2.5.  Pressure  history  for  a  radial  two-phase  shock  tube: 
Uniform  gas /liquid  mixture. 


Figure  3. 5. 2. 6  shows  the  spectral  analysis  of  the  limit-cycle  at  both  the 
centerline  and  the  upper  wall.  Both  the  MacCormack  scheme  and  the 
CFiAFT  upwind  procedure  predict  a  fundamental  frequency  of  approximately 
10.5  Hz.  However,  differences  exist  within  the  details  of  the  two  spectra. 
At  both  the  centerline  and  the  upper  wall,  the  upwind  solution  shows  a 
higher  percentage  of  energy  to  be  in  the  fundamental  tone,  while  the 
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MacCormack  solution  shows  more  energy  in  the  higher  overtones. 
Furthermore,  the  drop-off  in  the  magnitude  from  the  centerline  to  the  up¬ 
per  wall  is  much  larger  in  the  upwind  solution  when  compared  to  the  drop¬ 
off  in  the  spectra  of  the  MacCormack  solution.  Comparisons  of  Roe/TVD  and 
MacCormack  numerics  for  acoustic  boundary  layers  in  a  solid  propellant 
rocket  motor  also  show  comparable  differences  (see  Ref.  72). 
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3.5.3  Bt&lk  Liquid  Breakup 

Another  problem  of  fundamental  interest  which  can  be  studied  with 
the  gas-liquid  equilibrium  formulation  is  the  deformation  and  aerodynamic 
breakup  of  bulk  liquid.  Experimental  studies  have  recently  been  performed 
[73]  with  fly-out  of  bulk  liquid  to  provide  an  understanding  of  aerodynamic 
mechanisms.  The  liquid  fly-out  problem  is  of  present  interest  in  the  sce¬ 
nario  of  a  threat  missile  carrying  noxious  bulk  liquid  releasing  it  after  being 
"hit.".  Here  the  droplet  formation  mechanism  and  the  dispersion  of  the 
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droplet  cloud  is  of  critical  interest  to  the  defense  community.  In  the  exper¬ 
imental  studies,  a  cylindrical  pellet  of  liquid  flying  at  high  speeds  is  impul¬ 
sively  introduced  into,  stagnant  air.  The  problem  of  interest  is  how  the  bulk 
liquid  pellet  interacts  with  the  air  surrounding  it,  deforming  and  eventually 
having  droplets  being  stripped-off  from  the  interface  between  the  air  and 
the  liquid.  This  problem  directly  relates  to  instabilities  at  the  gas/liquid  in¬ 
terface  in  ETC/LP  and  LPG  problems  where  bulk  liquid  interacts  with  gas 
and  generates  droplets. 

The  schematic  of  the  problem  being  studied  is  shown  in  Fig.  3.5.3. 1.  A 
cylindrical  pellet  of  liquid  having  a  diameter  of  50  cm  and  a  length  of  50  cm 
is  impulsively  given  a  velocity  of  700  m/s.  The  flow  around  this  liquid  cylin¬ 
der  is  computed  using  a  dynamic  grid  which  translates  at  the  instantaneous 
velocity  of  the  grid  point  initially  located  at  the  liquid  center-of-mass.  The 
location  and  shape  of  the  liquid  mass  is  displayed  in  time  by  plotting  the 
liquid  mass  fraction  contours,  while  the  flowfield  around  and  within  the  liq¬ 
uid  mass  is  displayed  by  plotting  the  pressure  contours  in  the  flowfield. 


Fig.  3.5.3. 1.  Flowfield  initialization:  Flyout  problem. 


In  Figures  3.5.3.2a-3.5.3.2d  we  plot  the  liquid  mass  fraction  and  pres¬ 
sure  contours  at  the  following  four  times:  1,  2,  3,  and  4  ms.  Figure  3.5.3.2a 
shows  the  flowfield  at  1  ms.  At  1  ms  the  original  cylindrical  shape  (shown 
as  black  box)  begins  to  show  deformation  at  both  the  leading  edge  and  the 
trailing  edge.  At  the  leading  edge,  the  liquid  gets  compressed  while  in  the 
rear  the  liquid  at  the  comers  gets  stretched  out  in  the  form  of  fingers.  The 


pressure  contours  at  1  ms  show  the  classic  features  of  a  blunt  body  flow,  ex¬ 
cept  that  the  blunt  body  here  is  a  compressible  liquid  which  absorbs  part  of 
the  shock.  A  classical  bow  shock  is  formed  in  front  of  the  liquid,  and  a  wake 
with  a  neck  shock  is  formed  behind  the  liquid.  At  2  ms  (Fig.  3.5.3.2b),  the 
front  end  of  the  liquid  continues  to  get  compressed  and  pushed  out.  The 
rear  portion  of  the  liquid  shows  vortical  liquid  roll  up  with  the  gas.  The 
pressure  field  shows  that  the  bow  shock  continues  to  move  out,  while  the 
plume  expands  out  further  behind  the  liquid.  By  4ms,  the  liquid  compres¬ 
sion  results  in  two  large  fingers  of  liquid  being  pulled  out  at  the  top  and  bot¬ 
tom  surface  with  a  central  column  of  liquid  being  compressed  further.  The 
overall  deceleration  of  the  bulk  liquid  is  determined  in  the  course  of  the 
calculation  as  exhibited  in  Figure  3. 5. 3. 3.  The  liquid  slows  down  from  700 
m/s  to  250  m/s  in  the  4ms  of  this  prediction. 
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Fig.  3.5.3.3.  Velocity  history  of  bulk  liquid. 

The  above  calculations  highlight  the  ability  of  the  CRAFT  code  to  cap¬ 
ture  vortical  roll-ups  and  large-scale  turbulent  structures  which  produce 
mixing  between  phases  in  transient  multi-phase  flows.  To  take  the  problem 
of  gas-liquid  mixing  to  completion,  the  code  would  either  have  to  resolve 
very  fine  scales  (micron  level),  or,  would  have  to  be  coupled  to  a  droplet 
generation  model.  When  droplets  are  generated,  the  flowfield  would  consist 
of  a  gas-bulk  liquid  continuum  fluid  with  a  dispersed  phase  consisting  of 
droplets.  As  we  will  describe  in  Section  6,  our  work  in  solid  propellant  ETC 


(to  discretely  track  the  solid  propellant  particles)  has  yielded  methodology 
to  handle  such  liquid  flows  using  a  Lagrangian  solution  procedure  to  track 
the  droplets.  The  missing  link  in  developing  a  unified  code  to  solve  for  liq¬ 
uid  breakup  is  the  inclusion  of  a  dynamic  droplet  generation  correlation. 
Such  correlation  work  is  presently  under  development  in  the  threat  missile 
defense  community  and  should  have  significant  relevance  to  liquid  propel¬ 
lant  gun  simulation. 


4.0  SIMULATIONS  OF  ETC/LP  FIRING 


In  this  section,  we  describe  simulations  for  30  mm  liquid  propellant 
ETC  guns  using  the  gas-liquid  equilibrium  formulation  described  in  the  ear¬ 
lier  section.  We  begin  with  a  description  of  the  problem  set-up,  then  dis¬ 
cuss  the  simulation  of  Shot  17  fired  by  FMC.  In  the  final  section,  we  address 
unsteady  turbulence  issues  in  transient  flows  that  we  have  encountered  and 
dealt  with  in  our  ETC/LP  simulations. 

4.1  PROBLEM  SETUP 

The  gun  geometry  chosen  for  study  is  a  30  mm  ullage  tube  configura¬ 
tion  which  was  fired  at  FMC  and  is  denoted  as  "Shot  17."  Figure  4.1.1  shows 
a  schematic  of  the  gun  chamber  which  has  a  total  volume  of  167  cc  (length 
is  approximately  20.9  cm.  and  radius  is  approximately  1.59  cm).  Initially, 
the  plasma  occupies  the  center-core  ullage  tube  which  has  a  volume  of  19  cc 
and  extends  the  entire  length  of  the  chamber.  The  center-core  is  assumed 
to  burst  uniformly  along  its  entire  length  when  the  plasma  contained  within 
attains  a  critical  value  of  pressure  and  temperature.  The  burst  pressure  and 
temperature  of  the  plasma  are  specified  to  be  112  MPa  and  15067  K,  re¬ 
spectively.  These  conditions  were  obtained  from  calculations  performed  by 
Powell  using  his  one-dimensional  model  [70]. 
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Fig.  4.1.1.  FMC  Shot  17  ullage  tube  geometry  and  specifications. 


The  propellant  employed  for  the  Shot  17  firing  is  a  13.4  Molar  gelled 
solution  of  HAN/ Hydrocarbon.  The  total  mass  of  the  propellant  is  201  gm 
and  it  occupies  a  volume  of  148  cc.  As  we  described  in  Secbon  ,3.4,  a  two- 
step  combustion  model  is  utilized  in  CRAFT.  The  liquid  propellant  decom¬ 
poses  to  a  gaseous  state  and  then  combusts  to  form  the  final  products.  The 
details  of  this  two-step  combustion  process  and  the  respective  rates  for  the 
finite-rate  Arrhenius  form  are  given  in  Eqs.  (3.4.1)  and  (3.4.2).  The  numeri¬ 
cal  values  for  the  rate  coefficients  were  obtained  after  performing  one-di¬ 
mensional  numerical  experiments  to  determine  an  appropriate  rate  of  pres¬ 
sure  rise.  The  values  utilized  in  these  equations  are: 

Decomposition: 

A  =6.25x10^  B  =-4000.0 

V  V 

Combustion: 

A  =1.25x10^  B  =-4000.0 

c  c 

The  units  of  A  and  B  are  chosen  such  that  Ky  and  Kc  have  units  of  cc/(mole- 

sec). 


4.2  DISCUSSION  OF  SIMULATION 

The  results  from  our  simulation  of  the  30  mm  shot  are  compared  with 
the  experimental  measurements  obtained  by  FMC.  As  we  shall  discuss  in 
greater  detail  in  the  following  section,  the  transient  nature  of  the  flow  and 
the  mixing  generated  by  the  injection  of  high  pressure  plasma  are  expected 
to  generate  large  vortical,  turbulent  structures.  These  large  scale  structures 
are  captured  by  the  numerics  of  the  CRAFT  code.  Hence,  we  utilize  the  LES 
methodology  in  representing  the  turbulence  which  was  outlined  in  Section 
2,  wherein  only  the  small  scale  dissipative  turbulence  is  modelled.  For  the 
calculation  described  in  this  section,  an  algebraic  eddy-viscosity  SGS  model 
with  constant  length  scale  was  first  employed  to  model  the  effects  of  the 
small  scale  turbulence.  The  constant  length  scale  was  taken  to  be  1  percent 
of  the  radius  of  the  chamber.  Additional  sensitivity  studies  with  sub-grid 
turbulence  levels  and  other  related  issues  are  described  in  the  next  section. 
In  the  following  paragraphs,  we  discuss  and  compare  our  baseline  computa¬ 
tion  with  the  experimental  results. 
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The  measured  and  computed  pressure-time  history  on  the  chamber 
wall  at  position  A  (3.25  cm  from  breech)  is  shown  in  Fig.  4.2.1.  On  compar¬ 
ing  the  experimental .  data  with  the  numerical  pressure  curve,  we  observe 
that  the  initial  pressure  history  and  the  peak  pressure  value  is  reproduced 
closely  by  the  computations.  The  experimental  and  numerical  peak  pres¬ 
sure  value  is  approximately  290  MPa.  However,  the  experimental  peak  is  at¬ 
tained  at  a  time  of  1.5  ms  while  the  simulation  predicts  the  peak  to  be  at  2.0 
ms.  The  qualitative  nature  of  the  peak  is  represented  well  by  the  computa¬ 
tions.  The  data  shows  a  flat  pressure  peak;  the  peak  value  is  maintained  for 
about  0.5  ms.  The  computed  peak  also  has  a  flat  top  for  approximately  the 
same  time  period. 
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Fig.  4.2.1.  Pressure  history  at  Tap  A:  Experimental  data  vs  numerical 
computation  for  FMC  Shot  17. 

After  the  pressure  peaks,  the  velocity  of  the  projectile  continues  to  in¬ 
crease  until  it  equilibrates.  Correspondingly,  the  pressure  drops  to  its  equi¬ 
librium  value.  The  numerically  computed  equilibrium  pressure  level  is 
around  40  MPa  which  compares  well  with  the  experimental  equilibrium 
pressure  level.  However,  the  rate  of  drop  in  the  pressure  is  more  rapid  in 
the  computations  than  that  observed  experimentally.  The  reason  for  the 


faster  drop  in  computed  pressure  is  not  fully  understood.  One  possibility  is 
that  the  \drial  coefficients  of  the  gaseous  combustion  product  mixture  has  to 
be  better  represented  at  these  high  pressure  and  temperature  levels. 

The  pressure  histoiy  computed  at  Location  B  (17.5  cm  from  breech)  is 
shown  in  Figure  4.2.2  along  with  the  experimental  data.  The  general  obser¬ 
vations  made  earlier  for  Location  A  hold  here  as  well;  the  pressure  peak  val¬ 
ues  as  well  as  the  time  scales  match  well  with  the  data.  The  experimental 
peak  shown  here  was  obtained  at  an  azimuthal  angle  of  90  degrees.  This 
curve  differs  from  that  obtained  at  an  azimuthal  angle  of  0  degrees  which  has 
multiple  peaks.  The  experimental  data  reveals  significant  3-D  effects  which 
may  be  attributed  to  the  possibility  that  the  plasma  tube  did  not  rupture 
uniformly. 
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Fig.  4.2.2.  Pressiare  Mstoiy  at  Tap  B:  Experimental  data  vs  nnmerical 
computation  for  FMC  Shot  17. 

The  velocity  and  displacement  of  the  projectile  are  shown  in  Fig.  4.2.3. 
A  shot-start  pressure  of  20  MPa  was  specified  which  prevents  the  projectile 
from  moving  until  0.55  ms.  Once  the  projectile  overcomes  the  shot-start, 
the  velocity  becomes  proportional  to  the  integral  of  the  pressure  curve. 
Hence,  the  velocity  rises  exponentially  as  the  pressure  peaks  and  then  even- 


VEL  (IN  M/S) 


tually  equilibrates  to  a  value  of  approximately  1400  m/s.  Over  the  time  in¬ 
terval  of  the  computation,  the  projectile  travels  a  distance  of  approximately 
450  cm. 


(a)  VELOCITY  OF  PROJECTILE  (b)  DISPLACEMENT  OF  PROJECTILE 


TIME  (IN  MS)  (IN  MS) 

Fig.  4.2.3.  Numerical  predictions  for  projectile  velocity  (a) 
and  displacement  (b). 


4.3  UNSTEADY  TURBULENCE  ISSUES 

As  referred  to  earlier,  the  flowfield  within  the  ETC  flowfield  is  ex¬ 
pected  to  exhibit  large  scale  vortical  structures.  In  the  center  ullage  prob¬ 
lem,  vorticity  will  be  produced  by  the  projectile  and  enhanced  combustion 
should  occur  in  the  projectile  wake  region.  Consider  product  mass  fraction 
contours  in  the  chamber  (Fig.  4.3.1)  at  three  different  times,  which  are  good 
markers  for  mixing  effects  because  product  generation  increases  in  regions 
where  the  hot  plasma  mixes  effectively  with  the  liquid  propellant.  In  Fig. 
4.3.1,  the  blue  colors  at  the  top  and  bottom  of  the  chamber  denote  very  low 
product  concentrations  in  the  liquid  propellant.  As  we  proceed  towards  the 
core,  the  product  concentration  increases  sharply  across  the  interface  layer 
where  the  plasma/products  are  mixing  with  the  liquid.  Notice  the 
enhanced  mixing  in  the  projectile  wake  region.  To  highlight  the  enhanced 
generation  of  products  in  vortical  regions,  we  plot  instantaneous  contours  of 
mass  fraction  in  the  interval  1.395-1.495  (Fig.  4.3.2)  at  intervals  of  20  ps.  As 
Fig.  4.3.2  indicates,  the  product  generation  appears  to  predominantly  occur 
in  the  vortical  roll-ups  where  the  plasma  and  product  are  mixing  the  fastest. 
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Product  mass  fraction  contours  at  three  different  times  for  FMC  Shot  17  Simulation; 
1.75,  2.00,  and  2.25ms. 


Contours  of  combustion  product  mass  fraction  plotted  at  intervals  of  20u9 
beginning  at  1.395ms. 


The  large  vortical  structures  evidenced  in  Fig.  4.3.2  contain  the  large- 
scale  component  of  the  flow  turbulence.  To  illustrate  this,  we  process 
instantaneous  velocity  data  in  projectile  fixed  coordinates  for  the  time  pe¬ 
riod  0.95-1.4  ms.  From  the  instantaneous  velocity  values,  we  compute  the 
rms  value  of  the  fluctuating  velocity  component  at  various  axial  locations.  In 
Figs.  4.3.3  and  4,3.4  we  plot  the  radial  variation  of  the  rms  Reynolds  and 
Favre  averaged  fluctuations  at  positions  4  cm  and  17.5  cm  upstream  of  the 


Y  (IN  CMS)  time  (IN  MS) 


(a)  Radial  Variation  of  Instantaneous  Axial  Velocity  at  y  =  0.94  cms 

Figc  4,3.3.  Trabtaience  characteristics  in  tallage  tube  4.0cm  from  projectile 
base:  a)  radial  variation  of  Umisl  b)  instantaneous  axial  velocity 
at  y=0.94cm. 
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(a)  Radial  Variation  of  (b)  Instantaneous  Axial  Velocity  at  y  =  0.48  cms 

Fig.  4.3.4.  Turbulence  characteristics  in  ullage  tube  17.5cm  from  projectile 
base:  a)  radial  variation  of  Urmsl  instantaneous  axial  velocity 
at  y=0.48cm. 
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projectile.  The  Urms  peak  at  the  position  closer  to  the  projectile  (xp  - 
x=4cm)  is  closer  to  the  top  of  the  chamber  (y=0.9),  while  it  is  closer  to  the 
core  (y=0.4)  as  we  go  further  away  from  the  projectile  (xp  -  x  =  17.5  cm). 
These  observations  are  consistent  with  the  vortex  structures  evidenced  in 
Fig.  4.3.2.  Closer  to  the  projectile,  we  see  a  vortex  roll-up  extending  to  the 
top  of  the  chamber,  while  at  the  breech  end  the  vortex  structures  are  lim¬ 
ited  to  the  interface  shear  layer. 

Having  identified  that  large  scale  turbulent  structures  play  a  dominant 
role  in  the  combustion  process,  we  now  perform  numerical  experiments  to 
determine  the  sensitivity  of  the  pressure  curve  to  the  subgrid  turbulence 
levels  as  well  as  to  kinetic  rates.  We  repeated  this  calculation  with  rates 
which  are  twenty  percent  larger  than  the  baseline  rates  using  two  different 
subgrid  models.  Faster  rates  were  used  because  the  pressure  rise  in  the 
baseline  case  was  slower  than  observed  in  the  data.  The  first  subgrid  turbu¬ 
lence  model  used  was  the  same  as  the  baseline  case:  the  length  scale  is  one 
percent  of  the  radius  of  the  chamber.  The  second  subgrid  turbulence  model 
employed  is  the  Smagorinsky  formulation  with  Cg  =  .1.  We  designate  the 
baseline  case  as  Case  1  and  the  case  with  increased  rates  and  the  same  fixed 
length  scale,  as  Case  2.  Variants  of  Case  2  include  the  Smagorinsky  SGS 
model  (where  the  length  scale  expands  with  the  grid)  and  a  second  variant 
with  no  SGS  model. 

The  pressure  history  for  Case  2  is  plotted  in  Figures  4.3.5a  (x=3.25cm) 
and  4.3.5b  (x=17.5cm).  The  solid  line  corresponds  to  the  fixed  length  scale 
as  in  the  baseline  case,  while  the  dashed  line  corresponds  to  the 
Smagorinsky  subgrid  model  with  a  variable  grid/size  linked  length  scale. 
The  first  observabon  we  make  is  that  because  the  burning  rate  is  faster,  the 
pressure  rises  faster  than  the  baseline  case  and  reaches  a  value  of  290  MPa 
by  1.5  ms  as  in  the  experimental  data.  However,  having  reached  this  value, 
the  pressure  does  not  flatten  out  but  instead  rises  to  a  higher  level.  The 
second  observation  is  that  the  results  with  the  two  different  subgrid  turbu¬ 
lence  models  are  very  similar.  The  pressure  history  for  the  variable  length 
scale  (dashed  line)  flattens  out  a  little  more  than  the  fixed  length  scale 
(solid  line)  but  does  not  flatten  out  enough  to  prevent  the  pressure  from 
overshooting  the  peak  value  of  290  MPa.  The  shape  of  the  pressure  curve 
using  110  subgrid  scale  model  is  compared  to  that  with  a  subgrid  model  in 
Figure  4.3.6  (both  cases  have  the  same  burn  rates).  Without  the  subgrid 


model,  the  obsen^ed  "flattening"  does  not  occur.  These  results  are  a  prelim¬ 
inary  effort  to  assess  the  sensiti\dty  of  the  pressure  cur\^e  to  various  parame¬ 
ters.  Further  studies  on  grid  sensitivity  are  planned  and  we  intend  to  im¬ 
plement  more  advanced  LES  models  in  future  work. 
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Fig,  4,3, So  Sensitivity  of  pressure  history  to  bum  rates  and 
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5.0  CRAFT  CODE  EXTENSIONS  FOR  PRELIMINARY 
ETC/SP  SIMULATIONS 


In  this  section,  we  describe  simulations  of  solid-propellant  ETC  con¬ 
figurations  which  were  modelled  using  a  preliminary  fixed-bed  formulation. 
We  begin  by  describing  the  features  of  the  ETC/SP  designfollowed  by  a  dis¬ 
cussion  of  the  fixed  bed  formulation.  Simulations  for  two  shots  fired  by 
GDLS  are  discussed  and  compared  with  experimental  measurements. 

5.1  FEATURES  OF  THE  PROBLEM 

The  solid  propellant  ETC  gun  considered  in  our  study  had  a  center- 
core  ullage  tube  configuration  similar  to  the  ETC/LP  design.  A  schematic  of 
the  30  mm  ETC/SP  gun  analyzed  is  shown  in  Fig.  5.1.1.  As  in  the  ETC/LP 
case,  the  gun  chamber  has  an  ullage  tube  containing  plasma  which  bursts. 


Fig.  5.1.1.  Schematic  of  solid  propellant  ETC  gun  system. 

For  the  ETC/SP  case,  the  propellant  is  formed  of  spherical  balls  whose  ini¬ 
tial  diameter  is  approximately  1,000pm.  Loading  densities  are  high  and  the 
gas  porosity  is  of  the  order  of  0.4.  To  operate  with  these  high  loading  den¬ 
sities,  the  propellant  balls  are  deterred  chemically.  This  involves  diffusing  a 
chemical  deterrent  into  the  balls  such  that  burn  rate  increases  with  increas¬ 
ing  depth  from  the  surface  until  the  undeterred,  energetic  propellant  is 
reached.  A  typical  radial  variation  of  the  deterrent  profile,  which  is  shown 
in  Fig.  5.1.2,  indicates  that  the  properties  of  the  propellant  are  a  strong 
function  of  the  ball  diameter.  Simulation  of  the  solid  propellant  ETC  config¬ 
uration  requires  the  detailed  modelling  of  each  propellant  ball  as  well  as  the 
inter-phase,  non-equilibrium  processes  within  the  propellant  bed. 


Radial  Percsrit 

Fig.  5.1.2.  Representative  variation  of  deterrent  concentration  as  a  function 
of  deptihi  from  tine  surface  of  tine  ball  propellant, 

5.2  FIXED-BED  FOMMULATIOM 

As  a  first  step  in  the  implementation  of  the  two-phase,  nonequilibrium 
formulation,  preliminary  computations  were  performed  assuming  the  pro¬ 
pellant  bed  remains  fixed.  The  fixed  bed  assumption  eliminates  uncertainty 
associated  with  propellant  bed  deconsolidation  and  subsequent  fluidization. 
The  propellant  bed  boundary  is  instead  treated  as  a  regressing  internal 
boundary  through  which  combustion  products  enter  into  the  gas  core  (akin 
to  traditional  treatment  of  ablative  surfaces).  The  initial  volume  of  the  pro¬ 
pellant  bed  corresponds  to  the  solid  volume  of  the  propellant,  while  the 
ullage  volume  of  the  propellant  is  added  to  the  gas  core.  This  is  done  so  that 
the  plasma  being  pumped  in  has  additional  volume  for  pressure  relief,  and  it 
crudely  approximates  the  expansion  of  plasma  within  the  propellant  bed. 

The  numerical  methodology  of  the  fixed  bed  formulation  is  illustrated 
in  Tables  XII  and  XIII.  The  propellant  bed  is  divided  into  a  number  of  axial 
cells  and  the  diameter  of  propellant  balls  in  each  cell  is  tracked  in  time  as 
the  balls  combust  with  a  bum  rate  responding  to  the  pressure  in  the  gas 
core  at  the  propellant  bed  interface.  As  the  balls  burn  the  bed  regresses  and 
the  combustion  products  are  injected  into  the  gas  core.  The  regression  of 
the  bed  is  assumed  to  be  only  in  the  radial  direction  i.e.  the  axial  dimension 
of  the  propellant  cell  remains  unchanged  with  the  net  volume  change  being 
obtained  by  the  radial  movement  of  the  propellant  boundary. 
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Table  xn.  Coupling  of  Solid  Propellant  Bum  Model  Into  CRAFT 
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The  primary"  drawback  of  the  fixed  bed  formulation  is  its  inability  to 
integrate  through  the  propellant  bed.  Consequently,  the  physics  of  bed  de- 
consolidation  and  particle  entrainment  into  the  gas  core  cannot  be  mod¬ 
elled.  To  model  the  deconsolidation,  and  movement  of  the  propellant  balls 
due  to  interphase  drag,  a  more  general  equation  system  for  a  gas-particle 
mixture  is  required  where  the  particle  phase  can  occupy  significant  volume. 
Details  for  such  a  dense  two-phase  system  will  be  described  in  Section  6. 

5.3  SIMULATION  OF  GDLS  FIRINGS 

The  fixed  bed  formulation  was  used  to  simulate  GDLS  firings  of  the  30 
mm  ETC/SP  gun.  The  propellant  employed  for  these  firings  was  WC885.  In 
the  numerical  results  presented  here,  we  modeled  the  propellant  as  spheri¬ 
cal  balls  whose  diameter  is  specified  as  the  Sautered  mean  diameter  of  the 
propellant  balls  used  in  the  experimental  firings.  Thermodynamic  informa¬ 
tion  for  the  propellant  such  as  number  of  deterrent  layers,  chemical  energy 
release,  and  burn  rates  was  taken  from  closed  bomb /equilibrium  code  data 
provided  by  the  GDLS/OLIN  corporation.  We  note  that  the  deterrent  profile 
plays  a  crucial  role  in  determining  the  pressure  history  since  it  alters  both 
the  burn  rate  as  well  as  the  energy  release.  Figure  5.3.1  shows  a  representa¬ 
tive  ETC  calculation  performed  for  two  different  deterrent  layer  models  (4- 
layer  vs  15-layer),  but  otherwise  identical  conditions.  The  peak  pressure 
values  as  well  as  the  qualitative  nature  of  the  two  curves  show  significant 
differences  emphasizing  the  impact  of  the  deterrent  profile. 
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Fig.  5.3.1.  Sesasitiidty  of  pressnire  Mstoiy  to  deterreBt  model. 
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The  plasma  inflow  conditions  were  obtained  by  supplying  the  experi¬ 
mentally  measured  input  current  to  the  plasma  capillary  code  of  Powell  [70]. 
One  of  the  issues  not  fully  understood  was  the  plasma  burst  conditions  for 
the  moderator  tube.  The  burst  conditions  obtained  from  the  capillary  calcu¬ 
lation  appeared  too  high  for  the  material  of  the  moderator  tube.  It  was  de¬ 
cided  to  initialize  the  plasma  at  a  low  pressure  of  10  MPa  and  a  temperature 
of  15000  K.  The  reasoning  for  the  low  pressure  initialization  was  that  the 
plasma  might  burn  through  the  thin  moderator  tube  rather  than  burst 
through.  We  note  that  both  the  experimental  and  numerical  results  show 
sensitivity  to  the  plasma  initiation  processes  and  this  is  clearly  an  area  for 
further  study. 

Numerical  Simulations  were  carried  out  for  two  30  mm  shot  firings  by 
GDLS  —  Shot  45  and  Shot  122.  The  length  and  volume  of  the  gun  chamber 
are  approximately  20  cm  and  161  cc,  respectively,  and  the  mass  of  the  pro¬ 
jectile  is  328  gm.  The  mass  of  the  propellant  in  Shot  45  is  151  gm,  while 
Shot  122  had  a  lower  propellant  mass  of  116  gm.  The  current  inputs  to  the 
plasma  capillary  tube  for  the  two  shots  are  shown  in  Figs.  5.3.2  and  5.3.3  and 
indicate  that  Shot  45  has  significantly  more  plasma  energy  pumped  into 
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Fig*  5.3*2.  GDLS  Shot  45  input  data:  Current  history. 
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Fig.  5.3.3.  COLS  Shot  122  impiat  data:  Cmresit  history, 
the  gun  chamber.  This,  combined  with  the  higher  propellant  mass  in  Shot 
45  would  lead  us  to  expect  higher  performance  from  Shot  45  as  we  shall 
discuss  in  the  following  paragraph. 

The  numerically  computed  pressure  histoiy  and  the  experimental  data 
for  Shot  45  and  Shot  122  are  plotted  in  Fig.  5.3.4  and  Fig.  5.3.5,  respec- 
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NUMERICAL  SIMULATION 


Fig.  5.3.5.  GDLS  Shot  122  :  Comparison  of  computed  pressure  history 
with  experimental  data  at  Port  1. 

tively.  The  first  observation  we  make  is  that  the  peak  pressures  and  the 
risetimes  predicted  by  our  calculations  compare  well  with  the  experimental 
values.  After  the  pressure  peaks  the  pressure  drop  in  the  latter  half  of  the 
ballistic  cycle  is  not  predicted  as  well  by  the  computations.  This  is  probably 
because  the  projectile  is  moving  rapidly  at  this  point  and  the  fixed  bed  as¬ 
sumption  breaks  down.  The  good  agreement  between  the  computed  and 
experimental  data  also  extends  to  the  projectile  velocity  values.  The  exper¬ 
imental  muzzle  velocity  at  the  muzzle  exit  (130  cm]  for  Shot  122  is  857 
cm/s.  The  computed  muzzle  velocity  at  the  same  location  is  837  m/s  —  an 
error  of  2.4  percent. 

To  illustrate  the  complex  two-dimensional  flow  structure  within  the 
gun  chamber,  we  show  temperature  contours  at  three  different  times  in  Fig. 
5.3.6.  The  plasma  is  being  pumped  in  from  the  left  end  and  the  projectile  is 
moving  to  the  right.  At  the  earliest  time  level,  very  little  combustion  has  oc¬ 
curred  and  the  hot  plasma  fills  most  of  the  chamber.  As  time  increases,  the 
mass  of  the  products  increase  and  the  projectile  accelerates.  The  heavier 
products  push  the  lighter  plasma  out  and  a  complex  two-dimensional  struc¬ 
ture  develops.  In  the  last  snapshot,  the  projectile  velocity  is  considerable 
and  in  the  wake  behind  the  projectile,  periodic  vortex  structures  are 
evident. 
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Fig.  5.3.6.  Computed  log  temperature  contours  in  the  ETC/SP  gun  chamber  for  GDLS  Shot  45 
at  three  times:  0.65,  1.065,  and  1.465ms. 


6,0  MULTI-PHASE  FORMULATION  IN  CRAFT  FOR  GAS/PARTICLE 

SYSTEMS  IN  NON-EQUILIBRIUM 


In  this  section,  we  describe  a  generalized  formulation  for  a  mixture  of 
continuum  fluid  and  an  aggregate  of  incompressible  particles  with  the  two 
phases  being  out  of  equilibrium  with  each  other.  The  continuum  fluid,  in 
the  general  case,  can  itself  be  a  mixture  of  various  gas  species  and  bulk 
liquid  as  was  discussed  in  Section  3.  The  dispersed  phase  can  be  composed 
of  either  incompressible  solid  particulates  (as  in  solid  propellant  applica¬ 
tions),  or  liquid  droplets  (as  in  RLPG  applications).  The  generality  of  the 
approach  employed  for  dispersed  phase  physics  permits  utilizing  the  same 
non-equilibrium  framework  for  both  solid  and  liquid  propellant  problems. 
For  the  sake  of  clarity,  we  will  restrict  our  discussions  here  to  gas/solid 
particle  systems.  The  gas-phase  equations  in  CRAFT  are  solved  with  an 
Eulerian  procedure  using  the  upwind  numerical  framework  discussed  ear¬ 
lier.  Particle  motion  can  be  solved  using  either  an  Eulerian  or  a  Lagrangian 
procedure  with  each  procedure  entailing  its  own  set  of  advantages  and 
drawbacks.  Eulerian  methodology  has  been  implemented  in  CRAFT  for  solid 
propellant  and  rocket  plume/propulsive  flowfields  [10,11].  As  we  shall  dis¬ 
cuss  in  detail  later  in  this  section,  the  Lagrangian  procedure  is  more  appro¬ 
priate  for  the  ETC/SP  configuration.  We  present  numerous  validation  exer¬ 
cises  performed  with  the  new  Lagrangian  formulation  which  encompass  a 
wide  range  of  particle  volumetric  loadings  and  demonstrate  its  versatility. 
The  low  volumetric  loading  cases  have  applicability  to  rocket  nozzles,  mis¬ 
sile  plumes,  and  gun  muzzle  blasts,  while  the  high  volumetric  loadings  have 
applicability  to  gun  interior  ballistics,  fluidized  beds,  and  spray  applications. 

6.1  CONTINUUM-PHASE  EQUATIONS 

The  gas  phase  equations  for  dense  two-phase  flows  have  been  derived 
by  various  authors  for  applications  to  interior  ballisbcs,  sprays,  fluidized  bed 
flows,  etc.  The  equations  presented  here  are  analogous  to  those  originally 
developed  by  Cough  [74].  In  his  derivation.  Cough  has  averaged  the  equa¬ 
tions  representing  the  interaction  between  the  two  phases  over  a  region 
representing  a  cell  volume.  It  is  implicitly  assumed  that  cell  volumes  are 
larger  than  particle  volumes  in  this  formulation.  The  gas  phase  equations, 
cast  in  generalized  coordinates  are  given  below: 
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In  Eq.  (6.1.2),  pg  is  the  density  of  the  continuum  phase  which  is  taken  to  be 
a  pure  gas  to  simplify  our  notation.  We  note  that  this  is  merely  a  subset  of 
the  more  generalized  formulation  for  the  gas-bulk  liquid  mixture  described 
earlier.  If  the  continuum  phase  were  a  gas-liquid  mixture,  Pg  would  be  re¬ 
placed  by  pj^  where  p^^  is  defined  as  Pj^=Pg(l)g  +  PLitJ^-  The  term  Ug  in  Eq. 
(6.1.2),  denotes  the  volumetric  fracbon  of  the  continuum  phase  and  is  de¬ 
fined  as  ttg  =  1.0  -  Vp/V^gjj  where  Vp  is  the  total  volume  occupied  by  the  in¬ 
compressible  particle  phase  in  a  cell.  The  metric  quanbties, 
components  of  the  cell  face  normal,  L,  pointing  in  the  positive  ^  direction 
and  the  magnitude  of  the  metric  L  is  equal  to  the  cell  face  area.  For  multi¬ 
dimensional  flows,  the  vectors  F  and  G  have  forms  analogous  to  E  and  we 
similarly  define  metrics  M  and  N  corresponding  to  the  cell  faces  of  r|  and  C 
directions,  respectively.  The  quantity  U  denotes  the  volume  flux  through 
the  cell  face  in  the  ^  direction  and  is  defined  as  follows: 

U  =  + i.,v  + i.w  (6.1.3) 

X  y  z 


The  source  term,  S,  in  Eq.  (6.1.1)  contains  the  mass  transfer  terms 
due  to  combustion  of  the  solid  propellant  particles  (or  to  vaporization/ com¬ 
bustion  of  liquid  propellant  droplets)  and  is  given  as: 
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Here,  m^,  is  gas  generation  rate  from  particle  combustion  whose  functional 
form  will  be  defined  in  the  next  section.  The  vector,  contains  the 

nonequilibrium  drag  terms  and  is  given  as: 
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the  details  of  the  drag  coefficients  Ap  and  Bp  Eq.  (6.1.5)  have  been  taken 
from  correlations  for  packed  beds  (Gough,  Ref.  74). 

The  vector,  Fp,  contains  additional  source  terms  due  to  volumetric  ef¬ 
fects  and  is  given  as: 


In  Eq.  (6.1.6)  we  note  a  small  point  of  difference  in  our  methodology  with 
the  formulation  of  Gough  (Ref.  21).  Gough  represents  the  pressure  flux  term 
in  the  momentum  equation  as  a^VP.  For  numerical  convenience  in  retaining 
the  upwind  framework,  we  have  rearranged  this  term  as  follows: 

a  VP  =  v(a  p]-PVa  (6.1.7) 

In  Eq.  (6.1.7),  we  have  a  strongly  conservative  form  of  the  pressure  flux  term 
v(agP)  which  facilitates  the  implementation  of  an  upwind  numerical  proce¬ 
dure,  and,  we  discussed  earlier,  is  critical  to  capture  discontinuities  within 
the  flowfield.  We  note  that  the  equations  described  for  the  continuum  phase 
are  independent  of  the  formulation  utilized  for  the  particle -phase.  In  inter¬ 
acting  with  the  particle-phase,  the  gas  equations  require  only  the  inter¬ 
phase  interaction  terms,  namely:  particle  porosity  Up,  drag/heat  transfer 
terms  (F^),  and,  combustion  source  terms  (S).  Hence,  either  an  Eulerian  or 
a  Lagrangian  solution  procedure  for  the  particle  phase  can  be  coupled  as  a 
module  to  the  continuum  phase  equations  with  no  loss  of  generality  or 
physics. 

6.2  PARTICLE-PHASE  EQUATIONS 

In  formulating  the  numerical  procedure  to  be  implemented  for  the 
particulate  phase  in  ETTC/SP  simulations,  we  began  by  evaluating  the  relative 
merits  of  Eulerian  and  Lagrangian  formulations  for  the  particle  phase.  Our 
experiences  with  the  fixed-bed  formulation  (described  in  Section  5)  indi¬ 
cated  that  it  was  very  important  to  model  discrete  particulate  combustion 
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since  the  deterred  burn  properties  of  the  particle  are  a  strong  function  of 
the  particle  diameter.  The  diameter  of  individual  particles  can  be  accurately 
tracked  only  in  a  Lagrangian  scheme.  In  an  Eulerian  procedure,  scalcir  con¬ 
vection  equations  are  solved  for  the  surface  area  of  particles  in  each  particle 
group.  The  mean  diameter  in  each  group  would  then  have  to  be  coupled 
with  an  assumed  PDF  to  get  the  particle  diameter  distribution.  Hence,  an 
Eulerian  scheme  would  most  likely  not  provide  the  required  accuracy  in 
particle  combustion  problems  with  rates  sensitive  to  size  variation. 

Another  consideration  in  deciding  upon  the  numerical  scheme  is  the 
modeling  of  the  particle-parbcle  collisions.  For  dense  systems,  the  time  be¬ 
tween  collisions  decreases,  with  frequent  collisions  altering  the  local  volume 
fraction  occupied  by  particles.  This  directly  affects  the  gas-phase  equations 
since  it  alters  the  volume  available  for  the  gas  to  occupy.  In  a  Eulerian 
scheme,  while  the  interactions  between  different  particle  groups  can  be 
modeled,  the  particles  within  the  same  group  can  only  have  a  single  mean 
property  at  a  point  in  space.  This  condition  is  very  restrictive  when  there 
are  intersecting  streams  of  particles  because  the  Eulerian  scheme  will  yield 
a  single  "averaged"  stream  in  the  mean  direction  of  the  two  incident 
streams.  In  a  Lagrangian  scheme,  the  grazing  collisions  can  be  estimated 
with  stochastic  methods  although  the  methodology  is  quite  involved.  Based 
on  the  above  discussion,  a  Lagrangian  solution  procedure  was  chosen  as  the 
more  appropriate  methodology  for  ETC/SP  problems  of  interest  since  it  bet¬ 
ter  represents  the  discrete  physics,  and  can  be  extended  to  treat  collisions 
in  an  exact  manner. 

The  Lagrangian  equations  for  each  individual  particle  are  given  as  fol¬ 
lows: 

—  =  bp'^  (mass)  (6.2.1a) 

dt 

dQ  \  -V  -,/  \ 

m  — —  =  +A  [Q  -Q  )-v  VPh - a^ia  IVa  (momentum)  (6.2.1b) 

P  dt  P\  S  PI  P  N  „  P  \  SI  8  ' 
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m^-^  =  Bp^T^-Tpj  (energy)  (6.2. Ic) 

We  note  that  the  coefficients  for  the  pressure  based  burn  rates  and  the 
product  enthalpy  are  provided  from  experimental  data  at  various  points 


along  the  deterrent  profile.  Since  we  track  the  diameter  of  each  particle, 
the  coefficients  at  any  given  diameter  are  then  interpolated  from  the  deter¬ 
rent  data  curve.  A  second  observation  we  make  is  that  the  Lagrangian  equa¬ 
tions  described  above  do  not  contain  any  stochastic  modeling  for  collision 
effects.  This  is  an  area  for  future  study.  However,  we  have  included  an  in¬ 
tergranular  stress  term  for  densely  packed  beds.  From  an  Eulerian  view¬ 
point,  the  intergranular  stress  term  may  be  interpreted  as  a  macroscopic 
model  for  collisions.  The  form  of  this  intergranular  stress  has  been  adapted 
to  the  Lagrangian  formulation  from  the  term  given  by  Gough  [74]  for  his 
Eulerian  formulabon. 

We  again  note  that  an  Eulerian  particle  module  is  operational  in  CRAFT 
for  gas-particle  flows  where  the  particle  mass  loading  is  substantial,  but  vol¬ 
umetric  effects  are  negligible  [10,11].  Heavy  particles  such  as  AI2O3,  en¬ 
countered  in  rocket  propulsive  flows,  can  have  significant  mass  loadings 
(approximately  40%)  but  occupy  small  volumes.  The  particle  conservation 
equations  are  solved  in  an  upwind  fashion  to  mimic  the  gas-phase  numerics 
and  the  formulation  allows  for  non-equilibrium  between  phases  as  well  as 
phase-change  for  the  particles  as  they  heat-up  or  cool  down.  The  Eulerian 
package  has  been  used  extensively  to  study  transient  and  steady  flowfields  in 
rocket  nozzles,  exhaust  plumes,  and  other  aerodynamic  applications  where 
the  dilute  particle  formulation  is  appropriate. 


6.3  VALIDATION  OF  TWO-PHASE  FORMULATION  FOR  LOW 

VOLUMETRIC  LOADINGS 

The  two-phase  formulation  with  the  Lagrangian  solution  procedure  for 
the  particles  has  been  validated  by  computing  a  number  of  steady  and  tran¬ 
sient  unit  test  problems  for  low  volumetric  loading  and  comparing  these 
with  the  previously  developed  Eulerian  formulation  in  CRAFT.  The  steady 
state  problems  were  chosen  primarily  to  verify  that  particle  velocity  and 
temperature  equilibration,  as  well  as  non-equilibrium  interaction  with  the 
gas  were  implemented  properly.  The  transient  problems  were  computed  to 
verify  that  the  non-equilibrium  terms  alter  the  acoustic  propagation  in  the 
expected  fashion,  yielding  the  value  predicted  analytically  by  the  equilibrium 
two-phase  formulation  in  the  limit  when  non-equilibrium  effects  become 
negligible.  These  validation  cases  are  described  in  the  following  paragraphs. 
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6.3.1  Steady  Flow  Test  Cases 

The  steady-state  test  case  chosen  was  a  tube  with  1-D  gas-particle  flow 
within  it  (see  Case  I  in  Table  XIV).  The  gas  and  particle  flowfields  are  out  of 
equilibrium  at  the  inflow,  and  we  compare  the  velocity  and  temperature 
equilibration  lengths  as  they  travel  down  the  tube.  The  first  calculation  was 
performed  with  one-way  coupling,  i.e.,  only  the  particles  experience  the 
interphase  drag  while  the  gas  solution  remains  unaffected.  The  gas  velocity 
at  the  inlet  is  100  m/s  and  its  temperature  is  300  K.  The  particle  velocity 
and  temperature  are  half  that  of  the  gas  at  the  inlet.  The  particle  diameter 
is  lOgm  which  3delds  a  particle  lag  Reynolds  number  of  10.  The  particle 
volumetric  loading  is  negligibly  small  and  the  mass  loading  density  at  the 
inlet  is  0.5.  Since  the  volumetric  loading  is  small  the  Eulerian  particle 
package  for  dilute  two-phase  flows  is  used  to  generate  the  Eulerian  solution. 
The  solution  obtained  using  the  new  Lagrangian  package  is  compared  with 
the  Eulerian  formulation  which  has  been  validated  extensively  in  our  previ¬ 
ous  work  (Ref.  10). 

Table  XIV.  Validation  Studies  of  Lagrangian  Formulation 
for  Dense  Two-Phase  Flows 


CASEI:  STEADY  1-D  FLOW 

(LOW  VOLUMETRIC  LOADING) 

OBJECTIVE:  Verify  particulate  velocity 
and  temperature  equilibrium  with  gas 
flowfield  by  comparison  with  Eulerian 
solution 


CASE  II:  TRANSIENT  SHOCK  TUBE 

(LOW  VOLUMETRIC  LOADING) 

OBJECTIVE:  Verify  transient  interaction 
of  particle  motion  and  gas  flowfield  by 
comparison  with  Eulerian  solution _ 


CASE  III:  TRANSIENT  SHOCK  TUBE 

(HIGH  VOLUMETRIC  LOADING 
GUN  PRESSURE  CONDITION) 


OBJECTIVE:  Verify  transient  interaction 
of  particle  motion  and  gas  flowfield  by 
comparison  with  Eulerian  solution _ 
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In  Fig.  6.3. 1.1,  we  plot  the  velocity,  temperature  and  loading  density  of 
the  particles  as  they  travel  downstream.  In  Figs.  6.3.1.1a  and  6.3.1.1b 
(particle  velocity  and  temperature),  the  solid  line  is  the  Lagrangian  solution 
while  the  dotted  line  is  the  Eulerian  solution.  The  two  solutions  compare 
very  well  in  terms  of  the  equilibration  lengths  and  shape  of  the  curve.  The 
Lagrangian  solutions  for  the  velocity  and  temperature  are  smoother  than  the 
Eulerian  solution  because  these  quantities  are  integrated  continuously  in  the 
Lagrangian  case,  while  their  accuracy  is  limited  by  the  grid  size  in  the 
Eulerian  case. 
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- - Eulerian 


Tq  =  600 K 
Tp  >=3  300 K 


X  (IN  M) 


O.CCO  0.050  O.IDO  0.150  0.2C0 
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— ■ — — — --  Baseline  Injection  Rate 

„ — - 2x  Baseline  Rate 

- 4  X  Baseline  Rate 

*AII  density  levels  were  predicted 
to  be  the  same  (baseline  level). 

The  2x  and  4x  predictions  were 
shifted  for  clarity. 


Figisre  6.3. 1.1.  Steady-state  eqiriMbratioiffl  of  particles  (ome-way  coiipliiag). 

In  contrast  to  the  particle  velocity  and  temperature  profile,  the  parti¬ 
cle  loading  density  in  a  Lagrangian  formulation  shows  oscillations  (Fig. 
6.3.1.1c).  This  is  because  the  loading  density  in  the  Lagrangian  case  is  ob- 
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tained  after  the  time-step  integration  by  processing  the  positions  of  discrete 
particles  which  are  tracked  computationally.  Hence,  in  a  statistical  sense, 
we  are  approximating  a  continuous  solution  by  a  finite  sample  size. 
Therefore,  as  the  sample  size  becomes  larger,  we  expect  the  oscillations  to 
get  smaller.  This  effect  is  illustrated  in  Fig.  6.3.1.1c  which  shows  the  load¬ 
ing  density  computed  from  the  Lagrangian  calculation  for  three  different 
particle  injection  rates  at  the  inflow.  As  the  number  of  particles  being 
tracked  increases,  the  oscillations  which  were  observed  for  the  baseline  rate 
clearly  dampen  and  the  solution  becomes  smoother. 

The  1-D  gas-particle  flow  in  the  tube  is  now  computed  with  two  way 
coupling  between  the  gas  and  particle  phase,  i.e.,  the  gas  flowfield  is  affected 
by  the  interphase  drag  as  well.  The  gas  at  the  inflow  has  a  pressure  of  1  atm 
and  the  particle  mass  loading  is  0.5.  The  volumetric  loading  remains  negli¬ 
gibly  small.  The  initial  gas  velocity  is  400  m/s  while  the  particle  velocity  is 
half  that  at  200  m/s.  The  particle  lag  Reynolds  number  at  these  conditions 
is  65.  The  gas  and  particle  temperatures  are  identical  at  the  inflow. 
However,  this  changes  downstream  since  the  gas  heats  up  due  to  viscous 
dissipation  and  subsequently  this  results  in  interphase  heat  transfer. 

The  particle  and  gas  solutions  with  the  two-way  coupling  are  plotted 
in  Figs.  6.3. 1.2  and  6.3. 1.3  respectively.  These  solutions  are  plotted  2  ms  af¬ 
ter  particle  injection  began  at  the  inflow.  Prior  to  discussing  the  results,  we 
note  that  at  2  ms,  the  gas  solution  still  shows  transient  effects  because  the 
non-equilibrium  interactions  produce  waves  that  persist  after  equilibration. 
Figure  6.3. 1.2  shows  the  particle  solution  for  the  Lagrangian  and  Eulerian 
calculations.  The  particle  velocity  and  temperatures  compare  exactly  for  the 
two  cases.  The  particle  velocity  is  monotonically  equilibrating  with  the  gas 
velocity.  The  particle  temperature  rises  slightly  because,  as  mentioned  ear¬ 
lier,  the  gas  temperature  increases  due  to  viscous  dissipation.  The  particle 
loading  densities  at  2  ms  are  plotted  in  Fig.  6.3.1.2c.  The  particle  front  has 
propagated  to  a  distance  of  0.8  m.  The  Lagrangian  particle  front  preserves 
the  discontinuity  across  the  front,  while  the  Eulerian  solution  smears  the 
discontinuity.  The  Lagrangian  and  Eulerian  loading  density  profiles  are  dif¬ 
ferent,  as  expected,  because  the  Eulerian  procedure  treats  the  particle 
loading  density  as  a  continuous  function  while,  as  discussed  before,  the 
Lagrangian  solver  solves  for  the  individual  particles. 
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The  gas-phase  solution  at  2  ms  is  plotted  in  Fig.  6.3. 1.3.  We  observe 
that  although  the  particle  front  has  propagated  to  0.8  m  by  2  ms,  the  non¬ 
equilibrium  effects  on  the  gas  are  evident  only  to  a  distance  of  0.4  m. 
Beyond  0.4  m,  the  gas  velocity  and  pressure  are  the  original  initialized  val¬ 
ues.  In  the  region  where  the  gas  experiences  significant  drag,  the  gas  veloc¬ 
ity  drops  and  the  gas  pressure  rises.  In  fact,  the  gas  pressure  at  the  inlet 
has  gone  up  to  1.6  and  the  velocity  in  the  first  axial  cell  has  dropped  to  0.73 
from  the  initialized  value  of  1.0. 
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Figure  6.3. 1.3.  Gas  solution  2ms  into  injection  with  two-way  coupling. 
6.3.2  Transient  Test  Cases 

To  validate  the  Lagrangian  two-phase  code  for  transient  flowfields,  we 
set  up  a  two-phase  shock  tube  problem  (Case  II  and  III  in  Table  XIV).  The 
pressure  ratio  across  the  diaphragm  is  1 . 1 .  A  low  pressure  ratio  is  chosen  so 
that  we  can  estimate  the  acoustic  speed  rather  than  the  speed  associated 
with  a  strong  shock.  The  temperature  is  300°K  in  both  chambers.  The 
mass  loading  of  particles  is  0.5,  with  the  particle  diameter  being  10  mm. 
Two  sets  of  calculations  were  performed  with  different  pressure  levels.  The 
first  case  is  for  a  low  mean  pressure  of  1  atm.  At  this  pressure,  the  volu¬ 
metric  loading  of  particles  for  the  given  mass  loading  is  5  x  10*4.  The  sec¬ 
ond  case  is  for  a  high  pressure  of  117  MPa  which  is  in  the  regime  of  gun 
flowfields.  Since  the  gas  is  very  dense  at  this  high  pressure,  the  mass  load¬ 
ing  of  0.5  results  in  a  volumetric  loading  of  0.3. 
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Figures  6.3.2. 1  and  6. 3. 2. 2  show  the  gas  and  particle  solutions  for  the 
low  pressure  shock  tube  at  the  following  times:  0.4,  0.8,  1.2  ms.  The  gas 
pressure  profiles  are  plotted  in  Fig.  6.3.2.1a  and  the  gas  velocity  in  Fig. 
6.3.2.1b.  The  solid  line  shows  the  profiles  for  a  pure  gas  case  with  no  parti¬ 
cles  present.  The  two-phase  solutions  are  obtained  using  both  the  Eulerian 
and  the  Lagrangian  procedures  and  are  plotted  using  two  different  types  of 
dashed  lines.  However,  the  Lagrangian  and  Eulerian  solutions  fall  on  top  of 
each  other  and  appear  as  a  single  plot.  The  pressure  profiles  plotted  in  Fig. 
6.3.2.1a  indicate  that  at  these  low  volumetric  loadings,  the  acoustic  speed  of 
the  gas  slows  down  relative  to  the  pure  gas  phase  solution.  In  particular,  the 
pressure  profile  at  1.2  ms  indicates  that  in  addition  to  slowing  down,  the 
acoustic  wave  is  dispersing.  The  dispersive  nature  of  the  acoustic  waves  in 
non-equilibrium  flows  has  been  extensively  studied  in  the  literature  by  vari¬ 
ous  authors  [75]  for  linear  wave  equations  which  afford  analytical  solutions, 
and  is  a  well  known  phenomena.  The  decrease  in  the  acoustic  speed  at 
these  flow  conditions  can  also  be  verified  by  looking  at  the  equilibrium 
speed  of  sound  from  the  equilibrium  two-phase  formulation.  The  equilib¬ 
rium  acoustic  speed  when  the  second  medium  is  incompressible  is  as  fol¬ 
lows: 


where  is  the  acoustic  speed  for  the  mixture,  Cg  is  the  acoustic  speed  for 
a  pure  gas  phase  and  Gg  is  the  mass  loading  of  particles.  From  Eq.  (6.3.2. 1) 
it  is  clear  that  at  low  volumetric  loading,  (t)g  is  close  to  1.0  and  the  mixture 
acoustic  speed  drops  in  proportion  to  (1  + 

The  particle  loading  density  and  velocities  are  plotted  in  Fig.  6.3. 2. 2 
for  the  same  three  times.  The  particle  velocity  is  lower  than  the  gas  velocity 
(Fig.  6.3.2.1b)  which  shows  that  the  particles  have  not  yet  attained  equilib¬ 
rium.  Furthermore,  the  slope  of  the  particle  velocity  profile  on  either  side 
of  the  diaphragm  is  of  opposite  sign.  The  change  in  the  slope  of  the  particle 
velocity  around  the  original  diaphragm  location  is  reflected  in  the  particle 
loading  density  profiles  which  show  a  discontinuity  at  the  diaphragm  loca¬ 
tion.  On  the  right  side,  the  particle  loading  density  is  dropping  since  the 
gas  is  expanding,  while  on  the  left  side,  it  is  increasing  since  the  gas  is  be- 
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ing  compressed.  The  Lagrangian  and  the  Eulerian  loading  density  profiles 
compare  well  in  magnitude  and  shape.  However,  as  discussed  above,  the 
Lagrangian  solution  shows  oscillations  while  the  Eulerian  solution  is  smooth 
but  smeared.  We  note  that  since  the  parbcle  volumetric  loading  is  low,  the 
oscillations  in  the  parbcle  loading  density  are  not  reflected  in  the  gas  pres¬ 
sure  or  velocity  profiles.  However,  as  the  volumetric  loading  goes  up,  the 
oscillations  in  the  volumetric  loading  also  affect  the  gas  phase  solubon  as  we 
shall  describe  in  the  following  paragraphs  for  the  high-pressure  shock  tube. 
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Figure  6.3.2. 1.  Shock  tube  with  low  volumetric  loading:  Gas-phase  results. 
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The  solutions  for  the  high-pressure  shock  tube  (117  MPa)  are  show 
in  Fig.  6. 3. 2. 3  at  the  same  three  time  levels  as  in  the  earlier  case.  Figure 
6.3.2.3a  shows  the  pressure  profiles,  while  Fig.  6.3.2.3b  shows  the  corre¬ 
sponding  particle  loading  density  profiles.  The  first  observation  made  is 
that  the  acoustic  speed  goes  up  considerably  at  these  flow  conditions  and  is 
larger  than  the  pure-gas  phase  solution.  For  instance,  at  1.2  ms  the  pres¬ 
sure  wave  has  reached  the  end  walls  and  reflected  off,  while  the  pure  gas 
solution  in  Fig.  6.3.2.1a  has  stiU  not  yet  reached  the  end  walls.  The  increase 
in  acoustic  speed  can  also  be  deduced  from  the  equilibrium  sound  speed  ex¬ 
pression  in  Eq.  (6.3.2. 1).  Since  (j)g  drops  to  0.7,  the  volumetric  effect  of  (t)g^ 
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overrides  the  effect  of  the  mass  loading  and  increases  the  acoustic  speed. 
The  second  observation  made  is  in  regard  to  the  effect  of  particle  loading 
densities  on  the  pressure.  Since  the  volumetric  loading  is  substantial,  the 
oscillations  in  the  Lagrangian  loading  density  profiles  now  generate  ripples 
in  the  gas  pressure  profiles. 
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Figure  6.3.2.3.  Shock  tube  with  high  volumetric  loading:  Gas  pressure 
and  particle  loading  at  various  times. 
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6.4  CRAFT  VALIDATION  STUDIES  FOR  SOLID  PROPELLANT 

INTERIOR  BALLISTIC  APPLICATIONS 

In  this  section,  we  report  a  series  of  validation  exercises  which  were 
performed  to  evaluate  the  ability  of  CRAFT  to  model  solid-propellant  IB  ap¬ 
plications.  We  begin  by  computing  the  Love  and  Pidduck  Lagrange  problem 
to  ensure  that  the  gas-phase  thermodynamics  and  projectile  movement  is 
modelled  accurately.  Subsequently,  we  compute  a  series  of  cases  for  a  rep¬ 
resentative  30  mm  solid  propellant  gun.  These  calculations  are  done  in  se¬ 
quence  of  increasing  complexity  beginning  with  a  closed  bomb  case,  1-D 
calculations  for  different  propellant  loadings,  and  finally  culminating  with  a 
multi-dimensional  calculation,  with  plasma  injection,  which  is  representa¬ 
tive  of  a  ETC  gun.  For  the  closed  bomb  and  1-D  cases,  the  results  from 
CRAFT  are  compared  with  results  obtained  using  the  XKTC  code  [76]  which 
has  been  validated  extensively. 


6.4,1  Lagrange  Problem  (Love  and  Piddnck) 

In  Figure  6.4.1. 1  we  plot  the  pressure  distribution  at  various  times 
(0.4772ms-10.23ms)  for  the  problem  posed  by  Love  and  Pidduck  [77].  The 
pressure  profiles  compare  veiy  well  with  the  results  quoted  by  Love  and 
Pidduck.  The  differences  in  the  projectile  travel  and  base  pressure  for  each 
curve  is  shown  in  Table  XV.  The  maximum  difference  in  base  pressure  is 
0.43  percent  at  10.23  ms,  while  the  maximum  difference  for  projectile 
travel  is  -0.53  percent  at  5.154ms.  These  results  allow  us  conclude  that  the 
boundary  condition  at  the  accelerating  projectile  base  accurately  computes 
the  expansion  waves  emanating  from  the  projectile  without  introducing 
significant  conservation  errors.  We  note  that  Love  and  Pidduck  obtained 
their  results  using  a  method  of  characteristic  (MOC)  integration  procedure 
wherein  the  Riemann  variables  are  tracked  along  each  characteristic.  When 
there  are  waves  of  different  families  (expansion  and  compression)  crossing 
each  other,  the  accuracy  of  the  solution  would  depend  on  the  number  of 
characteristics  tracked.  Love  and  Pidduck  tracked  11  characteristics  in 
their  calculation.  It  is  interesting  to  observe  that  for  curve  4  (2.117  ms)  and 
curve  8  (7.137  ms)  the  pressure  error  flips  sign.  We  observe  that  at  these 
two  times,  the  reflection  from  the  breech  reaches  the  projectile  and  reflects 
back.  Using  a  larger  number  of  characteristics  may  have  an  impact  at  these 
times  when  a  MOC  procedure  is  employed. 
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Fig.  6.4. 1.1.  Pressure  profiles  at  various  times  computed  by  CRAFT  for 
Love  &  Pidduck  Lagrange  problem. 

Table  XV.  Love  &  Pidduck  Problem 


Projectile  Travel  (in  cm) 


Pressure  on  Projectile  (kg/cm^) 


Time 
(in  ms) 

Love 
&  Pidduck 

CRAFT 

Difference 

% 

Love 
&  Pidduck 

CRAFT 

Difference 

% 

0.4772 

2.4 

2.404 

-0.16% 

5651.3 

5652.6 

+0.023% 

0.9544 

9.28 

-0.02% 

5097.2 

5099.02 

+0.035% 

1.479 

21.4 

+0.317% 

4598.7 

4599.22 

+0.01% 

2.117 

42.191 

42.174 

4102.5 

4086.08 

-0.40% 

2.898 

75.4 

IHISI 

2970.3 

2975.30 

+0.16% 

3.859 

124.3 

-0.12% 

2161.6 

2165.32 

+0.17% 

5.154 

202.1 

-0.53% 

1535.2 

1541.06 

+0.38% 

7.137 

335.6 

335.545 

-0.016% 

1030.2 

1028.21 

-0.19% 

10.23 

571.9 

571 .84 

-0.0104% 

581.6 

584.104 

+0.43% 

6.4.2  Solid-Propellant  Interior  Ballistic  Calculations 

The  ability  of  CRAFT  to  accurately  model  conventional  solid-propellant 
guns  was  evaluated  by  simulating  a  30  mm  gun  for  various  propellant  load¬ 
ings.  The  geometry  of  the  gun  chamber  is  taken  to  be  cylindrical  with  a  vol¬ 
ume  of  141  cc.  The  propellant  employed  is  M-30  and  is  specified  to  be  in 
the  form  of  spherical  balls  with  a  diameter  of  1  mm.  The  mass  of  the  pro¬ 
jectile  is  specified  as  150  gm.  To  simplify  the  problem,  details  which  are 
secondary  for  the  purpose  of  code  validation  (i.e.  such  as  shot  start,  barrel 
resistance,  barrel  heat  loss,  etc.)  have  been  dropped  from  this  exercise. 

6.4o2ol  Closed  Bomb  Case 

The  pressure  and  temperature  history  for  the  closed  bomb  case 
(propellant  loading  0.25  g/cc)  are  shown  in  Fig.  6.4.2. 1.1.  The  CRAFT  and 
XKTC  solutions  match  exactly  on  the  scale  of  the  plot  (Figs.  6.4.2.1.1a  and  b). 
However,  when  the  scale  of  the  plot  is  blown  up  around  the  equilibrium 
value  minor  differences  can  be  observed  (as  seen  in  Figs.  6.4.2.1.1c  and  d). 
The  analytical  equilibrium  pressure  value  is  373.411  MPa  while  the  corre¬ 
sponding  temperature  value  is  2952.72  K.  The  CFIAFT  solution  reaches 
equilibrium  at  approximately  4.6  ms  with  a  pressure  of  373.374  MPa  and  a 
temperature  of  2952.43  K  which  are  in  excellent  agreement  with  the  analyt¬ 
ical  values.  Furthermore,  these  values  subsequently  remain  unchanged  indi¬ 
cating  that  the  code  is  conserving  energy  and  mass  exactly.  The  solution 
computed  by  XKTC  attains  a  equilibrium  value  of  373.374  MPa  and  2952.43 
K  by  4.6  ms.  However,  after  attaining  equilibrium  a  minor  amount  of  pres¬ 
sure  and  temperature  drop  is  observed. 


6.4.2.2  1-D  Gum  Calcialsitioms 

A  series  of  1-D  gun  tube  calculations  are  performed  for  propellant 
loadings  varying  from  0.1  g/cc  to  0.75  g/cc.  Here,  we  present  results  for 
two  cases:  the  0.1  and  0.75  g/cc  cases.  Fig.  6. 4. 2. 2.1  (a-d)  shows  the  breech 
pressure,  projectile  pressure,  projectile  velocity  and  displacement  respec¬ 
tively  for  the  CRAFT  (solid  line)  and  XKTC  (dashed  line)  solutions.  The  two 
solutions  compare  very  well  with  the  peak  breach  pressure  being  approxi¬ 
mately  10.5  MPa.  Furthermore,  the  projectile  velocity  and  displacement 
match  exactly  on  the  scale  of  the  plot.  The  corresponding  comparisons  for 
the  0.75  loading  case  are  shown  in  Fig.  6. 4. 2. 2. 2.  The  two  solutions  compare 
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1-D  gun  tube  calculation  for  0.75  g/cc  loading  case:  Pressure  history 
and  gun  performance  compared  for  CRAFT  and  XKTC  Solutions. 


remarkably  well  even  for  this  high  loading  case.  The  peak  breech  pressure 
is  around  370  MPa.  The  pressure  profile  for  the  XKTC  solution  shows  minor 
kinks  around  the  peak,  but  this  can  probably  be  eliminated  by  increased  grid 
resolution.  The  projectile  velocity  and  displacement  compare  very  well  with 
the  projectile  exit  velocity  being  approximately  1225  m/s. 

6.4.2.3  Multi-Dimensionail  ETC  Gim  Calculation  With  Plasma  Injection 

Having  validated  the  Lagrangian  formulation  in  CRAFT  for  1-D  calcu¬ 
lations,  we  simulated  a  representative  ETC/SP  configuration  to  demonstrate 
the  versatility  of  the  Lagangian  formulation  in  handling  complex  multi-di¬ 
mensional  flows  with  significant  spatial  variations.  The  gun  chamber  geome¬ 
try  is  identical  to  that  used  in  earlier  validation  studies  with  XKTC;  it  is 
cylindrical  with  a  diameter  of  30  mm  and  has  a  volume  of  141  cc.  The  aver¬ 
age  propellant  loading  in  the  chamber  is  0.75  g/cc.  However,  this  propel¬ 
lant  is  loaded  only  in  the  top  half  of  the  chamber,  and  therefore  the  local 
propellant  loading  is  above  1.0  g/cc.  The  central  portion  of  the  gun  cham¬ 
ber  contains  an  ullage  tube  into  which  plasma  is  pumped  in  eventually  caus¬ 
ing  the  tube  to  rupture  and  thereby  allowing  the  plasma  to  mix  with  the 
propellant  bed.  As  we  shall  describe  in  the  following  paragraph,  the 
Lagrangian  formulation  in  CFIAFT  allows  us  to  study  the  complex  process  of 
propellant  bed  deconsolidation  and  fluidization. 

Figures  6.4.2.3.1a  and  lb  show  the  porosity  and  log  temperature  con¬ 
tours  at  0.36  ms.  The  temperature  plot  shows  that  the  plasma,  at  this  time, 
is  generally  limited  to  the  central  core  of  the  chamber  except  at  the  pro¬ 
jectile  end  where  it  reflects  back  in  the  form  of  a  spherical  wave.  The  re¬ 
flection  of  the  plasma  from  the  projectile  end  is  reflected  in  the  porosity 
contours  of  the  propellant  bed  which  indicate  that  the  bed  deforms  to  ac¬ 
commodate  the  plasma.  The  propellant  bed  also  shows  a  slight  bulge  at  the 
inflow  end  where  the  high  pressure  plasma  is  trying  to  expand  out  against 
the  propellant  bed.  The  corresponding  temperature  and  porosity  contours 
at  0.61  and  0.86  ms  are  shown  in  Figs.  6. 4. 2. 3. 2  and  6. 4. 2. 3. 3  respectively. 
At  0.61  ms,  the  deconsolidation  and  fluidization  of  the  propellant  bed  be¬ 
come  apparent.  The  propellant  descends  towards  the  center  of  the  cham¬ 
ber  and  tries  to  follow  the  projectile  which  is  accelerating  forward.  The 
non-equilibrium  between  the  two  phases  is  highlighted  by  the  fact  that  the 
propellant  lags  the  projectile.  The  log  temperature  contours  show  that  the 


temperature  in  most  of  the  chamber  has  dropped  considerably  since  sub¬ 
stantial  amounts  of  cool  (relative  to  plasma  temperatures)  products  are  being 
produced  by  the  burning  propellant.  Furthermore,  the  two-dimensional 
structure  of  the  temperature  contours  indicate  that  some  high  temperature 
plasma  is  still  present  in  the  core  as  well  as  the  projectile  end.  In  fact  at 
the  projectile  end,  the  heavy  products  being  formed  in  the  chamber  push 
the  much  lighter  plasma  up  all  the  way  to  the  top  wall.  These  effects  con¬ 
tinue  to  be  evident  at  0.86  ms  (Fig.  6. 4. 2. 3. 3),  at  which  point  the  projectile 
has  traveled  further  ahead  with  most  of  the  propellant  already  burnt. 
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loading  case  at  0.36ms;  Loading  density  contours. 
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7.0  FUTURE  DIRECTION 


The  ETC  gun  simulation  modeling  described  in  this  report  has  pro¬ 
vided  an  enhanced  level  of  understanding  of  fundamental  processes  in  liquid 
and  solid  propellant  interior  ballistics.  The  sophisticated  upwind /implicit 
computational  framework  utilized  in  CRAFT  has  provided  the  ability  to  cap¬ 
ture  strong  discontinuities  in  a  non-oscillatory  manner,  and,  to  treat  acous¬ 
tic  disturbances  with  minimal  numerical  attenuation.  The  up\vind/ implicit 
framework  has  been  systematically  extended  to  include  advanced  thermo¬ 
chemical  features  (virial  EOS,  liquid  EOS,  gas/liquid  equilibrated  formula¬ 
tion)  and  Eulerian  and  Lagrangian  dispersed-phase  solutions  applicable  to 
both  liquid  droplet  and  solid  particle  analyses.  The  gas-phase  equations 
have  been  extended  to  account  for  particle /droplet  volumetric  effects  and 
the  Lagrangian.  formulation  was  utilized  to  simulate  packed-bed  solid  propel¬ 
lant  interior  ballistic  flows  with  great  success. 

In  its  present  form,  CFIAFT  has  all  the  necessary  ingredients  for  the 
generalized  analysis  of  liquid  and  solid  propellant  interior  ballistic  flows.  In 
addition  to  the  ETC  problems  described  in  this  report,  ram  accelerator 
(RAMAC)  [25-28]  and  LPG  [7]  problems  have  been  analyzed.  Specialized  grid 
methodology  (d5niamic  gridding,  patching/blanking  and  embedding  proce¬ 
dures)  and  the  utilization  of  LES  procedures  for  representing  the  transient 
turbulent  eddy  structure  are  common  to  all  the  interior  ballistic  flowfields 
analyzed.  Current  work  on  extending  CRAFT  to  utilize  hybrid  struc¬ 
tured/unstructured  grid  procedures  [78],  and,  on  the  improvement  of  sub- 
grid  stress  models  for  combusting  flows  [79]  will  greatly  enhance  its  basic 
framework. 

For  interior  ballistic  flows,  future  work  should  focus  on  thermochemi¬ 
cal  details  and  on  coupling  with  thermal  and  structural  solvers  to  obtain  ac¬ 
curate  transient  boundary  conditions  incorporating  wall  temperature  varia¬ 
tions  and  vibrations.  The  thermochemical  work  must  focus  on  first-princi¬ 
ples  methodology  to  replace  heuristic  modeling.  For  liquid  propellant 
problems,  the  details  of  what  actually  occurs  at  the  gas/liquid  interface  has 
not  been  modeled.  There  is  a  need  to  focus  on  the  droplet  formation  prob¬ 
lem  from  a  first-principles  viewpoint  since  the  mechanism  relates  directly 
to  dynamic  aspects  of  Icirge  eddy  motion  at  the  interfaee,  which  is  being  di¬ 
rectly  simulated  using  LES  methodology.  In  proposed  work  for  ERDEC  re- 
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lated  to  bulk  liquid  flyout/breakup,  a  first-principles  investigation  of  this 
problem  will  be  initiated  which  is  of  direct  relevance  to  liquid  propellant  IB 
flows . 

For  solid  propellant  IB  flows,  the  Lagrangian  framework  provides  gen¬ 
erality  but  is  numerically  complex.  Automated  procedures  to  control  the 
number  of  particles  that  must  be  tracked,  and,  improved  numerics  to  deal 
with  smoothing  cell-to-cell  transitions  of  particles  and  in-cell  averaging  pro¬ 
cedures  are  required.  In  addition,  a  first-principles  approach  to  deal  with 
particle/particle  collisions  is  needed  to  eliminate  ad  hoc  "fluidized-bed"  ex¬ 
tensions  to  conventional  drag/heat  transfer  laws. 
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APPENDIX  A 


INVISCID  FLUX  JACOBIAN  FOR  GAS-LIQUID  EQUILIBRIUM  FORMULATION* 
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where  f  represent  the  component  of  the  cell  face  area  in  the  ^  direction. 

U  is  the  contravariant  velocity  and  is  defmed  as 
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The  derivatives  appearing  in  the  matrix  are  defined  as  follows: 
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APPENDIX  B 


EIGENVALUES  AND  EIGENVECTORS  FOR  THE  GAS-LIQUID  FORMULATION* 


The  eigenvalues  of  the  Jacobian  A  are  represented  as 
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here,  Cj  and  Cl  are  the  isothermal  sound  speeds  for  the  i**^  gas  species,  and  the  liquid  species,  respectively. 


The  left  and  right  eigenvectors  are  represented  as  follows 
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*  Generalized  L  and  R  matrix  formulation  from  Ref.  3. 
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In  the  above  eigenvectors  i  =  i^j  +  is  the  normal  of  the  cell  face  in  the  |  direction.  The  vectors  m  andw 

are  two  arbitrary,  mutually  perpendicular  vectors  to  the  vector  i  .  The  terms  U  ^  F,  W ,  are  the  dot  product  of  the 
velocity  vector  with  these  three  unit  vectors. 
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J  KENNEDY 
MN38-3300 

10400  YELLOW  CIRCLE  DR 
MINNETONKA  MN  55343 

2  OLIN  ORDNANCE 

ATTN  V  MCDONALD  LIBRARY 
HUGH  MCELROY 
PO  BOX  222 
ST  MARKS  FL  32355 

1  PAUL  GOUGH  ASSOCIATES  INC 

ATTN  P  S  GOUGH 
1048  SOUTH  ST 
PORTSMOUTH  NH  0380L5423 

1  PHYSICS  INTERNATIONAL  COMPANY 
ATTN  LIBRARY  H  W  WAMPLER 
2700  MERCED  ST 

SAN  LEANDRO  CA  94577 

2  ROCKWELL  INTERNATIONAL 
ATTN  BA08  J  E  FLANAGAN 
JGRAY 

ROCKETDYNE  DIVISION 
6622  CANOGA  AVE 
CANOGA  PARK  CA  9134 

1  PRINCETON  COMBUSTION  RESEARCH  LAB 
ATTN  N  MESSINA 

PRINCETON  CORPORATE  PLAZA 
1 1  DEERPARK  DRIVE 
BLDG  IV  SUITE  119 
MONMOUTH  JUNCTION  NJ  08852 

2  SCIENCE  APPLICATIONS  INC 
ATTN  J  BATTEH 

L  THORNHILL 

1519  JOHNSON  FERRY  RD 

SUITE  300 

MARIETTA  GA  30062-6438 


1  ELI  FREEDMAN  &  ASSOCIATES 

ATTN  E  FREEDMAN 
2411  DIANA  RD 
BALTIMORE  MD  21209 

1  VERITAY  TECHNOLOGY  INC 

4845  MILLERSPORT  HWY 
PO  BOX  305 

EAST  AMHERST  NY  14051-0305 

1  BATTELLE 

TWSTIAC 
505  KING  AVE 
COLUMBUS  OH  43201-2693 

1  GENERAL  ELECTRIC  CO 

ATTN  DR  J  MANDZY 
DEFENSE  SYSTEMS  DIVISION 
MAIL  DROP  43  220 
100  PLASTICS  AVE 
PITTSFIELD  MA  01201 

1  STATE  UNIVERSITY  OF  NEW  YORK 

DR  W  J  SARGEANT 
DEPT  OF  ELECTRICAL  ENGNRNG 
BONNER  HALL  ROOM  312 
BUFFALO  NY  14260 

1  OFHCE  OF  NAVAL  RESEARCH 

ATTN  DR  S  G  LEKOUDIS 
800  NORTH  QUINCY  STREET 
ARLINGTON  VA  22217 

1  US  ARMY  MISSILE  COMMAND 

ATTN  AMSMI  RD  SS  AT  DR  B  J  WALKER 
REDSTONE  ARSENAL  AL  35898 

1  WRIGHT  LABORATORY 

ATTN  WL  POPS  DR  BALU  SAKAR 
WRIGHT  PAT  AFB  OH  45433 

1  PRATT  &  WHITNEY 

ATTN  DR  S  SYED 
MAIL  STOP  715  89 
PO  BOX  meoo 

W  PALM  BEACH  FL  33410-96(X) 

1  SDI 

DR  L  H  CAVENY 

INNOVATIVE  SCIENCE  AND  TECH 
THE  PENTAGON 
WASHINGTON  DC  20301-71(X) 
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1  UNITED  DEFENSE  LP 

ATTN  DR  J  H  KOO 
ARMAMENT  SYSTEMS  DIV 
MS  M443 

4800  EAST  RIVER  ROAD 
MINNEAPOLIS  MN  55421 

1  US  ARMY  ERDEC 

ATTN  SCBRD  RTM  DR  M  MILLER 
APG  MD  21010-5423 

1  DEPARTMENT  OF  THE  AIR  FORCE 

ATTN  AFOSR  DR  J  M  TISHKOFF 
BLDG  410  ROOM  237A  C  WING 
BOLLING  AFB  DC  20332 

1  NAVAL  SURFACE  WARFARE  CENTER 

ATTN  DR  G  MOORE 
DAHLGREN  DIVISION 
DAHLGREN  VA  22448 

1  NASA  LANGLEY  RESEARCH  CENTER 

ATTN  DR  J  M  SEINER 
MAIL  STOP  165 
HAMPTON  VA  23681-0001 

1  PHILLIPS  LAB  ASTRONAUTICS  LAB 

ATTN  OLAC  PL  RKFT  DR  J  LEVINE 
BUILDING  8350 
EDWARDS  AFB  CA  93523 

1  INSTITUTE  FOR  APPLIED  TECHNOLOGY 

ATTN  DR  W  G  REINECKE 
4030  2  WEST  BREAKER  LANE 
AUSTIN  TX  78759 

1  US  ARMY  CORP  OF  ENGINEERS 

ATTN  DR  J  P  BALSARA 
WATERWAYS  EXPERIMENT  STATION 
3909  HALLS  FERRY  ROAD 
VICKSBURG,  MS  39180-6199 

1  DEPUTY  COMMANDER 

ATTN  CSSD  WD  L  DR  R  BECKER 
US  ARMY  SP  AND  STRTGC  DEFNS  CMND 
PO  BOX  1500 
HUNTSVILLE  AL  35807 


1  US  ARMY  RESEARCH  OFFICE 

ATTN  DR  D  MANN 
ENGINEERING  SCIENCES  DIVISION 
4300  S  MIAMI  BLVD 
RSRCH  TRI  PK  NC  27709 

1  WRIGHT  LABORATORY 

ATTN  WL  FIMG  DR  J  SHANG 
WRIGHT  PAT  AFB  OH  45433 

1  INSTITUTE  FOR  DEFENSE  ANALYSIS 

ATTN  DR  H  WOLFHARD 
1801  N  BEAUREGARD  STREET 
ALEXANDRIA  VA  22311 

1  PHILLIPS  LAB  ASTRONAUTICS  LAB 

ATTN  OLAC  PL  RKFT  DR  P  KESSEL 
BUILDING  8350 
EDWARDS  AFB  CA  93523 

1  GENERAL  SCIENCES  INC 

ATTN  DR  P  ZAVITSANOS 
205  SCHOOLHOUSE  ROAD 
SOUDERTON  PA  18964 

1  NASA  MARSHALL  SPACE  AND  FLIGHT  CENTER 

ATTN  DR  P  MCCONNAUGHEY 
MAIL  CODE  ED  32 
MSFC  AL  35812 

1  NAVAL  AIR  WARFARE  CENTER 

ATTN  DR  K  C  SCHADOW 
WEAPONS  DIVISION 
CODE  C02392 

PROPULSION  RESEARCH  BRANCH 
CHINA  LAKE  CA  93555-6001 

1  NAVAL  AIR  WARFARE  CENTER 

ATTN  DR  K  A  GREEN 
AIRCRAFT  DIVISION 
OFFICE  OF  SCIENCE  AND  TECHNOLOGY 
CODE  4.0  MS  70 
WARMISTER  PA  18974-0591 

1  DEPUTY  DIRECTOR  AEROPROPULSION 

ATTN  H  M  BRILLIANT 
WRIGHT  LABORATORY 
WRIGHT  PAT  AFB  OH  45433 
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1 
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ORGANIZATION 

PETC 

ATTN  DR  D  WELDMAN 
COAL  COMBUSTION  R&D 
PITTSBURGH  PA  15236 

ARMY  HIGH  PERFORM  CMPTG  RSRCH  CTR 
ATTN  DR  T  TEZDUYAR 
UNIVERSITY  OF  MINNESOTA 
1 100  WASHINGTON  AVE  S 
SUITE  101 

MINNEAPOLIS  MN  55415 

NASA  AMES  RESEARCH  CENTER 

ATTN  DR  D  COOPER 

NAS 

MS  258  5 

MOFFET  FIELD  CA  94035 

TRW  DEFENSE  SYSTEMS  GROUP 

ATTN  SR  T  C  LIN 

BALLISTIC  MISSILE  DIVISION 

PO  BOX  1310 

BLDG  527  RM  706 

SAN  BERNARDINO,  CA  92402 

NAVAL  POST  GRADUATE  SCHOOL 
ATTN  DR  D  NETZGER 
DEPT  OF  AERONAUTICS 
HALLIGAN  HALL 
BUILDING  234 
MONTEREY  CA  93943 

PETC 

ATTN  DR  J  EKMAN 

DIRECTOR  COAL  COMBUSTION  R&D 

PITTSBURGH  PA  15236 

COMBUSTION  RESEARCH  &  FLOW 

TECHNOLOGY  INC 

174N  MAIN  ST 

BLDG  3  PO  BOX  1150 

DUBLIN  PA  18917 


NO.  OF 

COPIES  ORGANIZATION 

ABERDEEN  PROVING  GROUND 

4  CDR  USACSTA 

ATTN:  S  WALTON 
G  RICE 
D  LACEY 
CHERUD 

1  DIR  USAHEL 

ATTN  J  WEISZ 

24  DIR  USARL 

ATTN  AMSRL-WT-PA, 

T  MINOR 
G  WREN  (2  CP) 

W  OBERLE  (2  CP) 
PTRAN 
J  DESPDUTO 
T  COFFEE 
MNUSCA 
G  KELLER 
D  KOOKER 
F  ROBBINS 
R  ANDERSON 
G  KATULKA 
K  WHITE 
P  CONROY 
S  RAY 
A  JUHASZ 
AMSRL-WT-PB, 

E  SCHMIDT 
P  PLOSTINS 

AMSRL-WT-PC,  R  FTFER 
AMSRL-WT-PD,  B  BURNS 
AMSRL-WT-WD,  J  POWELL 
AMSRL-WT-WG,  P  KASTE 


DR  W  J  SARGEANT 
DEPT  OF  ELCTRCL  ENGNRNG 
BONNER  HALL  RM  312 
BUFFALO  NY  14260 
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2  RARDE 

ATTN  DR  C  WOODLEY 

DR  G  COOK 

GS2  DIVISION 

BUILDING  R31 

FORT  HALSTEAD 

SEVENOAKS  KENT  TN14  7BP 

ENGLAND 

1  WEAPONS  SYSTEM  DIVISION 

ATTN  DR  ANNA  WILDEGGER-GAISSMAIER 
PO  BOX  1500 

SALISBURY,  SOUTH  AUSTRALIA  5108 
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Intentionally  left  blank. 


DIST-8 


USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  pubUshes.  Your  comments/answers 
to  the  items/questions  below  will  aid  us  in  our  efforts. 

1.  ARL  Report  Number  ARL-CR-240 _ Date  of  Report  August  1995 _ 

2.  Date  Report  Received _ ^ _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  puipose,  related  project,  or  other  area  of  interest  for  which  the  report 

will  be  used.)  _ _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of  ideas,  etc.) 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  doUars  saved,  operating  costs 
avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to 
organization,  technical  content,  format,  etc.) _ 


Organization 


CURRENT  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the 
Old  or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 
(DO  NOT  STAPLE) 


OFRCIAL  BUSINESS 


WO  POSTAGE 
NECESSARY 
IF  MAILED 

IN  THE 

UNITED  STATES 


